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ABSTRACT 


This work accomplishes signal denoising using the Bootstrap method when the 
additive noise is Gaussian. The noisy signal is separated into frequency bands using the 
Fourier or Wavelet transform. Each frequency band is tested for Gaussianity by 
evaluating the kurtosis. The Bootstrap method is used to increase the reliability of the 
kurtosis estimate. Noise effects are minimized using a hard or soft thresholding scheme 
on the frequency bands that were estimated to be Gaussian. The recovered signal is 
obtained by applying the appropriate inverse transform to the modified frequency bands 
The denoising scheme is tested using three test signals. Results show that EFT-based 
denoising schemes perform better than WT-based denoising schemes on the stationary 
sinusoidal signals, whereas WT-based schemes outperform FFT-based schemes on 
chirp type signals. Results also show that hard thresholding never outperforms soft 
thresholding; at best its performance is similar to soft thresholding. 
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EXECUTIVE SUMMARY 


The main goal of this thesis was to develop a denoising scheme to recover signals 
that are distorted by additive white Gaussian noise, regardless of the signal’s sperctral 
content. Two approaches were investigated, a short-time Fourier transform and a wavelet 
transform-based approach. Both transforms are used to decompose a signal into several 
frequency bands and to min im ize the noise in each band. Each frequency band content is 
tested for Gaussianity, which is accomphshed by investigating the signal’s kurtosis. It is 
well known that the signal kurtosis is equal to three only when the signal is Gaussian with 
zero mean. Therefore this parameter can help to provide a rehable determination whether 
or not a data sequence is Gaussian. However, the kurtosis estimate of a short data 
segment is not always rehable. The Bootstrap method is used in the kurtosis estimation to 
overcome this difficulty. The Bootstrap is a statistical scheme, which improves the 
rehabihty of a parameter estimate in situations where conventional techniques are not 
valid because of short data length issues. The Bootstrap uses data sequences that have the 
same length as the original signal, but are obtained by randomly resampling the data 
using replacement. These resampling is done many times and each data set is treated as 
repeated experiment. Then, the parameter of interest is estimated for each of these 
resampled sequences, to obtain a statistic for the parameter of interest. This resulting 
statistic can then be used to perform hypothesis tests on the parameter value. 

If the band is estimated to be Gaussian (i.e., the noise is dominant), then 
thresholding is apphed to the band to min im ize the noise effects. Two thresholding 
techniques are considered, hard and soft threshold. Hard thresholding coefficients are 0 
for Gaussian data and 1 for non-Gaussian data. Whereas soft thresholding coefficients are 
obtained by considering the data’s closeness to Gaussianity. The closer to being Gaussian 
the band specific data is, the smaller the soft thresholding coefficients are and vice versa. 
The denoised signal is obtained by applying the appropriate inverse transformation. 

Three different test signal types are selected to investigate the performances: a 
sinusoid, a chirp with constant amphtude, and a chirp with an RC time constant- tik e 


XV 



amplitude increase. The mean square error (MSB) and a distance measure defined 
between original and recovered signals are selected to compare performances. Results 
show FFT-based denoising schemes perform better than WT-based denoising schemes on 
the stationary sinusoid signal type, whereas WT-based schemes outperform FFT-based 
schemes on chirp type signals. Finally, results show that the soft thresholding scheme 
always performs at least as well as or better than the hard thresholding one. 



1. INTRODUCTION 


In many data transmission and storage systems, noise gets introduced into data, 
which reduces the signal quality. As a result, various filtering techniques have been 
designed to suppress noise contributions in order to improve the overall signal quality. 

Fourier and Wavelet transforms decompose a noisy signal into several frequency 
bands. Traditional filter design methods have requirements on the frequency, magnitude 
and phase of the signal such as passband ripple, stopband attenuation, transition width, 
and phase constraints. The assumption behind these design criteria is that the signal is 
restricted to be in a certain frequency band and that the frequencies outside this specific 
band are treated as distortion. Note that this paradigm breaks down when signal and 
distortion terms overlap in frequency. 

Denoising attempts to remove the noise and to recover the original signal 
regardless of the signal’s frequency content. The basic idea is to look at each frequency 
band of interest and to minimize its noise effect by retaining the dominant component. 

The band is left untouched when the signal is dominant so as not to lose the signal 
component, while thresholding is apphed when the noise is dominant. 

This thesis discusses a denoising scheme that implements frequency band specific 
thresholding schemes using the Bootstrap method and the kurtosis. Chapter 2 presents the 
processing techniques and the Bootstrap method. Chapter 3 discusses the proposed 
denoising scheme. Simulation results are presented in Chapter 4. Finally, Chapter 5 
presents conclusions. 


1 
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n. BACKGROUND 


A. PROCESSING TECHNIQUES 

Signal processing allows users to extract relevant information from a given signal. 

When the raw data does not allow extracting the desired information, transformations to 
another domain may be performed to do so. The most common transformation types are 
the Fourier and Wavelet transforms, which are discussed in the following two sections. 

I. Fourier Analysis 

Fourier analysis allows the representation of a given signal as a hnear 
combination of complex sinusoids with different frequencies. For periodic signals this is 
called the Fourier Series. This representation becomes extended to the Fourier 
Transform when the signal is aperiodic. Both representations are discussed in the 
following subsections. 

a Fourier Series 

A periodic signal x{t) with period may be represented as an infinite 
hnear combination of complex exponentials: 

’ ( 2 . 1 ) 

k=-oo 

271 

jk—t 

where /g = l/Tg is fundamental frequency, e is caUed the harmonic, and at is 

the weight. This representation is called the Fourier series representation [1], where 
the set of coefficients are cahed the Fourier series coefficients which are obtained by 

1 r -jk—t 

a^ = — ^x{t)e dt, k = . (2.2) 

2o 7 -„ 

The coefficients a/t’s are a measure of the strength of the signal’s components at the k^^ 
harmonic of the fundamental frequency. 
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b. Fourier Transform 

For aperiodic signals the Fourier series representation can be extended. 

Assuming that (t) is a periodic signal with a period , and x(t) represents one period 

of Xp (t), then as increases x^ (t) is identical to x(t) over a longer interval. Therefore, 
in the limiting case [2] 

x(t) = lim x„ (0 • (2.3) 

Replacing the lim its of the integral in equation (2.2) by - 0 ° and 0 °, and multiplying both 
sides by , leads to the Fourier Transform of x(t) given by 

X{f)=]x{t)e-^^^^dt, (2.4) 


where f = klT^ and X{f) = Tpa^. The Inverse Fourier Transform is used to recover the 
original signal x(t) from X(f), and is defined as 

x(t) = ]x(f)e^^^^‘df. (2.5) 


signals as 


The Discrete Time Fourier Transform (DTFT) is defined for discrete-time 


X((o) — ^x(n)e , 

n=oo 


( 2 . 6 ) 


and the corresponding inverse is given by 


x(n) =— [ X (co) e^"d(0 . 
2k •’ 


2n 


(2.7) 


It should be noted that the DTFT is periodic with period 27t . 

The Discrete Fourier Transform (DFT) for a finite time discrete-time 
signal x(n) with n=0,N-I is defined as 

X(yt) = ^x(n)e-^'"""^. (2.8) 

n=0 
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Its inverse is given by 


(2.9) 

The Fast Fourier Transform is a computationally efficient implementation of the 
Discrete Fourier Transform when the signal length is a power of two. The Discrete 
Fourier Transform requires N multiphcations whereas the Fast Fourier Transform 
requires (A^/2) log 2 N multiplications. 

c. Short-Time Fourier Transform 

The basic Fourier Transform allows the passage from the time domain to 
the frequency domain. However, one may need to preserve both the time and the 
frequency information contained in non-stationary signals. Unfortunately, the basic 
Fourier Transform does not provide such dual information. 

Time loca liz ation can be introduced by windowing the signal before using 
the Fourier Transform, where the window size is selected short enough to assume the 
signal inside the window is stationary. Then, the Fourier transform may be implemented 
on each windowed signal portion. The resulting transformation is called the Short-Time 
Fourier Transform and expressed as 

5(x, /) = j x{t)w(t-x) e-^^^^'dt , (2.10) 

t 

where w{t-x ) denotes the sliding window centered around x . Note that ^(x , /) 
provides a two-dimensional representation of the signal frequency information at various 
times X and is a “local” specfrum of the signal x{t) around the analysis point x [3]. 

Many different window types may be selected, depending on the characteristics of the 
Short-Time Fourier Transform desired. In addition, the window length determines the 
resolution in time and in frequency; good time resolution requires a short window, 
whereas good frequency resolution requires a long window. The joint time-frequency 
resolution of the Short-Time Fourier Transform is lim ited by the uncertainty principle. A 
short time window results in a loss of frequency resolution, and vice versa. 
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2. Wavelet Analysis 

A wavelet is defined as a small wave which has its energy concentrated in time 
and frequency [4]. Such a characteristic is useful for the analysis of transient, non¬ 
stationary or time-varying phenomena. By comparison, sinusoidal functions used in 
Fourier Analysis have a constant amphtude. 

The Wavelet Transform (WT) provides an alternative to the Short-Time Fourier 
Transform (S lFl ), as it uses short windows at high frequencies and long windows at low 
frequencies [5], while the S lFl uses windows of constant size, as illustrated in Figure 1. 
Note that the Wavelet Transform still satisfies the uncertainty principle, however, the 
time resolution becomes arbitrarily good at high frequencies, while the frequency 
resolution becomes arbitrarily good at low frequencies. The following subsections 
discuss different types of Wavelet Transform that are used for continuous or discrete time 
signals. 


frequency frequency 




(a) STFT (b) WT 


Figure 1. (a) Time-Frequency resolution for the STFT, (b) Time-Frequency resolution for the WT. 
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a. The Continuous Time Wavelet Transform 

The Continuous Time Wavelet Transform may be applied to continuous 
signals and two dimensions in the transfer domain are continuous. It is similar to the 
Fourier Transform in that, it is obtained by projecting the signal onto a basis function. 
However, the Continuous Time Wavelet Transform projects the signal onto scaled and 
shifted versions of the wavelet function while the Fourier Transform uses complex 
exponentials as basis functions. The Continuous Time Wavelet Transform is defined as 

C(5,x) = ^ ? x(t)\^dt, 

where \|/ (t) is the wavelet function, x denotes the shift in time and s is the scale factor 
that denotes compression or expansion in time. The inverse transform for finite K is 
obtained by 

x(t) = ^[ ■\c(s,x )t|/ (-——)dsdx , 
where the parameter K is given by 



and 'F (to) is the Fourier Transform of the wavelet function v|/ (t). [6] 
b. The Discrete Time Wavelet Transform 

The Discrete Time Wavelet Transform is the discrete time version of the 
Continuous Wavelet Transform. It is used for discrete time signals and the dimensions of 
the transform domain are discrete as well. The Discrete Wavelet Transform is obtained as 

C(a ,b) = V-^x(n)\|/ , 

n fa a 

where a, b and n are the discrete parameter versions of s, X and t given in Equation 
(2.11), respectively. The scahng factor a has another restriction as a = where 


( 2 . 11 ) 


( 2 . 12 ) 


(2.13) 


(2.14) 
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j=0,l,. . .,\ogi{signal length). The common choice for is 2, as it allows for fast 
algorithms. 

c. Mallafs Algorithm 

The Discrete Wavelet Transform, can be implemented by using Mahat’s 
algorithm when 0 ^ = 2 [4], and is illustrated in Figure 2 for a three-level decomposition. 
The signal is passed through a high-pass and a low-pass filter, both of which have a 
bandwidth of half the signal spectmm. Then, following the Nyquisf s mle, the outputs of 
the filters are subsampled by two. The subsampled high-pass filter output is called the 
Detail Sequence and the subsampled low-pass filter output is called the Approximation 
Sequence. This procedure may be recursively applied to the approximation sequence 
obtained at previous levels. Note that a signal of length 2^ can be decomposed only j 
times, because the approximation sequence has only one sample left after j levels. 



Figure 2. Mallat’s algorithm. 


B. HIGHER ORDER STATISTICS 
I. Moments and Cumulants 

The first four moments for a real valued and stationary signal x{n) are given by 


m^ = E[x{n )^, 

(2.15) 

ty^) = E{x(n)xin+x^)], 

(2.16) 


8 







,X 2 ) = E\x{n)x{n +x^)x{n+x 2 )], 


(2.17) 


m4(Xi ,X2,X3) = £'{x(n)x(n+Xj)x(n+X2)x(n+X3)] . (2.18) 

The first two moments are equal to the mean and correlation functions respectively. The 
first four cumulants are given by 

Ci = mi, (2.19) 

C2(Xi)=m2^i)-mj(Xj)\ (2.20) 

C3(Xi ,X2) = m3(Xi ,X2) -m^irn^ (xj + -Xj)] +2(m^)^ (2.21) 

^4 (^ 1 5^2 3 ) ^ 4 ^^ 1 ’^3 ) m2(Xj )w^ (X 3 X 2) tn2(X 2)^^ (^3 ^ l) 

-m2(X3)«^(X2 -Xi)-mj[m5(X2-Xi,X3 -Xj) 

+ m3(Xi,X2) + m3(Xi,X3) + m3(X2,X3)] (2.22) 

+ (m,)^[m2(Xj) + m2(X2) + m2(X3) + m2(X3 -X^) 

+ m 2 (X 3 -X 2 ) + m 2 (X 2 -Xj)] - 6 (mj)‘*. 

Note that the second and third order cumulants are identical to the second and third order 
moments respectively when x(n) is a zero-mean process. 

Cumulants have properties that make them more desirable than moments. Some 
of these properties are 

a. Each cumulant is independent of all lower order cumulants. 

b. AH cumulants of order greater than two are equal to zero for Gaussian 
processes. Hence, any Gaussian process is completely characterized by its first two 
cumulants. Therefore higher-order cumulants can be used to estimate the degree of non- 
Gaussianity of a process. 

c. Cumulants of the sum of two independent statistical processes are equal to 
the sum of their respective cumulants. [6] 

2. Variance, Skewness and Kurtosis Measures 

Setting Xj,X 2 ,X 3 equal to 0 in (2.20), (2.21), (2.22), and assuming that m7=0, leads 

to the variance , skewness Y 2 , and kurtosis y^ measures: 

9 



y^^E\x{nf] =C 2 ( 0 ), 


(2.23) 


(2.24) 

(2.25) 

where a is the standard deviation around the mean of the signal. [7] 

C. THE BOOTSTRAP 

The Bootstrap is a powerful technique for assessing the accuracy of a parameter 
estimator in situations where conventional techniques are not applicable [8]. In many 
apphcations, one needs to estimate one or more parameters of a random process, and/or 
calculate some statistical parameters such as the mean or variance. Most of the estimation 
techniques used for this purpose assume that the set of samples used in the estimation is 
large enough to reach asymptotic results. However, in practice this assumption usually 
does not hold as the sample set may not be large enough or the samples may be non¬ 
stationary. The Bootstrap scheme randomly reassigns the observations, recomputes the 
estimates many times, and treats these reassignments as repeated experiments. In this 
section, first the basic Bootstrap principle is stated. Next, the usage of Bootstrap in 
estimating the confidence interval for a parameter is discussed, and finally this discussion 
is extended to hypothesis testing. 

1. Basic Principle 

Let X2,..., Xnj be a collection of n independent and identically distributed 

random variables drawn from an unknown distribution, D. Let p denote an unknown 

statistical parameter of D such as the mean or the variance, and p denote an estimator of 
p, calculated from x. If the estimate p is to be used in place of the real parameter p, it 

may be important to know the samphng distribution of p . The distribution of p may be 
estimated with the Bootstrap method, and is obtained by resamphng many times from a 

distribution D , chosen to be close to D, such that D approaches D as n ^ oo. Note that 
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the choice of D is not unique. If the type of D is known but its statistical parameters are 

not known, then D is chosen as a distribution of the same type as D, with the statistical 
parameters obtained from x. For example, if we know that the data is Gaussian but do not 
know its mean and variance, we perform the resamplings assuming the data to have the 
same mean and variance as x has. This approach is called the parametric Bootstrap. If 
nothing is known about D, the resamphngs are drawn from x with replacement, so that 
each value in a resample set has probability equal to 1/n. This approach is called the 
nonparametric Bootstrap. 

2. Parameter Confidence Interval 

The Bootstrap principle may also be applied to obtain a (l-a)100% confidence 
interval for the parameter p. First, the data set is resampled many times such that each 
resampled set is of the same size as the data. Next, an estimator is obtained from each 
resampled set, where k =\,...,N and N is the number of repetitions. Then, the estimates 

are sorted in increasing order. Finally the indices of the lower lim it p^ and the upper 
limit pjj of the estimates in 

P(p^<p<p^) = l-a 


are obtained by using 


L = 


Na 


and 


U =N-L+\, 


where [aJ denotes the integerpart of the value A.[8] 

3. Hypothesis Testing 
a. Description 

Suppose one needs to perform a hypothesis test, such as H : p< p^^ 
against the hypothesis K: p> p^, where p^ is given. A new statistic, defined as 


(2.26) 


(2.27) 


(2.28) 
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p-p. 


(2.29) 


may be selected, where o is the estimator of the standard deviation a of p . The 
estimator for the standard deviation o of p wiU be defined later. 

To perform the hypothesis testing, one first draws random sequences 
Xj*, X 2 , of the same size as x, with replacement from x. Note that * does not denote 

complex conjugation, but means that this is a resampled set, or a parameter obtained from 

* 

a resampled set. Then the statistic T is estimated from each sequence x as 


P -P 


(2.30) 


where p and a are estimated parameters obtained from the resample x*, instead of x, 
and the constant p^ is replaced with p . Note that o is included as a scale factor in the 

calculation of f . Dividing by a is called Bootstrap pivoting and it is done to ensure f 
is asymptotically pivotal when n ^ 00 , i.e., the asymptotic distribution of f does not 
depend on any unknown parameters. Replacing p^ with p and using Bootstrap pivoting 

is important because the Bootstrap distribution of T -(p - p) I a is a better 
approximation to the distribution of T = (p- Pq)/o under H, than the Bootstrap 
distribution of S* - p -p is to the distribution of S =p -p^ under H [8, 9]. Next, the 
set of test statistics T\,T 2 ,...,T n are sorted by increasing order, and the hypothesis H is 


rejected if T >TiM), where M is chosen according to N and the level of significance a as 


[8, pp. 62] 


M = (N+\)(\-a). 


(2.31) 
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Note that the test statistie T is given by 


P-Po 

a 


(2.32) 


when the hypotheses to be tested we, H : p = against K : p^ p^. 

b. Estimation of the Standard Deviation for p 

The parameter a ean be estimated by using the Bootstrap. Towards that 
end, resamples of the same length as the data set v, called x* are drawn from x randomly 
with replacement, to obtain a total of B resamples. After resampling, the Bootstrap 

estimates p are calculated in the same manner as p was, but with resamples x* instead 
of X. As a result, the standard deviation o of p is estimated by 

In this chapter we presented the processing techniques that we used in the 
proposed denoising scheme. In the next chapter we will discuss the proposed denoising 
scheme. 
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m. DENOISING USING THE KURTOSIS AND THE BOOTSTRAP 


A. INTRODUCTION 

Fourier and Wavelet transforms decompose a signal into several frequency bands. 
As a result, white Gaussian noise components affect all frequency bands. The proposed 
technique examines each frequency band and tends to minimize the white noise 
contribution. First, we obtain the signal time-frequency representation (or time-scale 
representation in the Wavelet transform case). Next, each frequency band is tested for 
Gaussianity, and thresholding is performed on the spectral location found to be Gaussian. 
Last, the signal is transformed to the time domain using the appropriate inverse 
transform. The following sections discuss two different implementations using Fourier 
and Wavelet transformations. 

B. FAST FOURIER TRANSFORM BASED DENOISING 

The FFT-based denoising scheme is performed in four steps, as illustrated in 
Figure 3. For computational convenience, the signal length is assumed to be 512 points. 

A longer data set could have been used, resulting in more frequency bands, and more 
computations. A data set of minimum length of 512 allows separation in time and 
frequency. A shorter segment will not provide a reasonable number of time-frequency 
cells. 



Figure 3. FFT-based denoising scheme. 
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1. Short-Time Fourier Transform Step 

In the first step, the data is divided into 32 data-point segments and weighted by a 
triangular window. A four-to-one overlap is used to obtain a reasonable data size. Then, 
each segment is transformed to the frequency domain with the Fast Fourier Transform. 

This results in 17 frequency bands. A larger data set would have allowed for a larger 
segment size. The STFT information is contained in a matrix of dimensions 61 by 32 
elements for a data length of 512 points, where the first number corresponds to the time 
dimension and the second one corresponds to the frequency dimension. 

2. Gaussianity Test Step 

The second step tests the Gaussianity of the transformed values in each frequency 
band, where real and imaginary parts of the data are tested separately with the kurtosis. 
Recall that the norma liz ed kurtosis value for a Gaussian data sequence is equal to 3. 

Hence, looking at the kurtosis value may give an idea of a sequence’s Gaussianity, 
provided the data length is sufficient for the estimation to be meaningful. However, the 
sequences in this particular scheme are of length equal to 61, which may not be sufficient 
for meaningful estimations. Therefore, the Gaussianity test is implemented as a 
hypothesis test, with the hypothesis H -.kurtosisof thedata = 3 , and the hypothesis 
K -.kurtosisof thedata ^ 3 . The sequence is found to be Gaussian when the hypothesis H 

is hue, meaning the kurtosis value for the data sequence is 3 within the specified 
confidence interval. However, the sequence is found to be non-Gaussian when hypothesis 
K is true, meaning the kurtosis is not equal to 3 within the selected confidence interval. 

This hypothesis testing is performed using the Bootstrap method discussed earher in 
Chapter 2. An empirical confidence interval a =0.05 was selected here. 

3. Thresholding Step 

Signal components are found in a particular sequence when some frequency 
bands are estimated to be non-Gaussian for both real and imaginary parts. Thresholding is 
apphed to each band estimated to be Gaussian, to minimize noise effects in these bands. 
Two thresholding schemes are considered below. 
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a. 


Hard Thresholding 


Hard thresholding zeroes out ah values in the frequency band that is found 
to be Gaussian. The hard thresholding coefficient is 

f 0, if thebandisGaussian 

Ch=\ 

1, if thebandisnotGaussian. 


(3.1) 


b. Soft Thresholding 


Soft Thresholding is obtained by multiplying values in the specific 
frequency band that is found to be Gaussian, by a coefficient between 0 and 1. Using a 
coefficient of 0 is the same as hard thresholding, whereas using a coefficient of 1 is the 
same as leaving the frequency band undisturbed. The soft thresholding coefficient is 
calculated using 



(3.2) 


where Y 4 is the bootstrapped kurtosis of the particular frequency band, and | | denotes 

the absolute value. The bootstrapped kurtosis value Y 4 is l im ited not to exceed 4.5, which 

will be explained later. The bootstrapped kurtosis value is obtained by using the 
Bootstrap principle, which calls for resamplings from the data set many times with 
replacement, to obtain N resampled sets of the same length as the original data set. Next, 

the kurtosis value for each resample is found. Finally, the Bootstrapped kurtosis Y 4 is 
defined as the estimated mean obtained from kurtosis values. It should be noted that the 
thresholding coefficient c is a function of the frequency band’s degree of Gaussianity. 
Equation (3.1) shows that the coefficient c gets closer to 0 as the bootstrapped kurtosis 

value Y 4 for a specific frequency band gets closer to the theoretical value 3, and vice 
versa. Therefore, the closer a frequency band gets to being Gaussian, the smaller 
contribution it has after soft thresholding. 

Recall that frequency bands to be thresholded are those that passed the 
Gaussianity test. Therefore, one would expect their corresponding bootstrapped kurtosis 


17 



estimates to be close to the theoretical value 3. However in some cases, the estimated 
bootstrapped kurtosis value may be far off from 3. For example, a few values in the band 
may be higher than the others especially when a frequency bin contains a short data 
segment. In this case, the band may pass the Gaussianity test as most of the band is 
Gaussian, but may still have a bootstrapped kurtosis value high enough to obtain a 
thresholding coefficient c greater than 1. Note that a thresholding coefficient greater than 
1 would amphfy noise contributions in the frequency band, instead of suppressing them. 

Therefore, the bootstrapped kurtosis value is lim ited to 4.5 in (3.1) to avoid such 
potential noise amphfication. This value is chosen empirically. Thus, if the bootstrapped 
kurtosis value of a specific frequency band is greater than 4.5, the data point with the 
largest absolute value in that band is stored away and replaced with a zero. Then, the 
bootstrapped kurtosis for that band is estimated again. This procedure is repeated until a 
kurtosis value less than or equal to 4.5 is obtained. However, the number of repetitions 
should be l im ited to ensure that most of the data is not zeroed out, to prevent 
bootstrapped kurtosis estimation problems. An empirical lim it of one third of the number 
of time frames defined in the STFT. Note that if the estimated kurtosis value is still 
greater than 4.5 after that many iterations, then it is set to 4.5, which provides a 
thresholding coefficient equal to 1. This causes no change on the frequency band under 
consideration. 

After multiplying the frequency band with the thresholding coefficient, the 
values that were removed to lim it the kurtosis estimation are reinserted into their original 
locations. This allows for the non-Gaussian values, which may correspond to the signal 
components, to not be affected by the thresholding step. 

4. Inverse Fourier Transform 

Real and imaginary parts of the frequency bands are combined to form the time- 
frequency representation matrix after the thresholding step, and the recovered signal is 
obtained using the inverse Fourier transform. 
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C. WAVELET TRANSFORM BASED DENOISING 

The Wavelet transform-based denoising scheme is similar to the FFT-based 
scheme, as illustrated in Figure 4. Note that only real valued transforms are obtained in 
this case due to the real wavelet transform kernel form selected in our work (Daubechies 
wavelets of order 3,4 or 5). The choice of Daubechies wavelets was motivated by the 
earher works in GSM signal denoising reported by Aktas and Mantis [10, 11]. This 
denoising scheme is performed in four steps. 



Figure 4. WT-based denoising scheme. 


1. The Wavelet Transform 

First, the Wavelet transform is apphed to the noisy signal. Three-, four-, or five- 
level decompositions are considered in this work, because simulations showed that 
higher-level decompositions do not improve performance, given the data length 
considered (512). 

2. Gaussianity Test 

The approximation and detail coefficients of a Gaussian data set remain Gaussian 
[12]. Thus, detail and approximation coefficients are tested for Gaussianity and a decision 
is made for each. 
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3. Thresholding 

The same thresholding procedure as that was considered earlier for the FFT-based 
denoising scheme is apphed to the detail coefficients. However, thresholding the 
approximation coefficients may result in removing a significant portion of the signal. One 
can use one of the following four schemes to threshold the approximation coefficients to 
minimize potential distortion: Apply the thresholding scheme as that used on the detail 
coefficients, leave coefficients unperturbed, filter the approximation coefficients using a 
median filter of order three, or use a predictor of order two. The last two options were 
investigated as ways to smooth the approximation coefficients but they did not perform 
as good as the first two. Results for these operations are discussed in the next chapter. 

4. Inverse Wavelet Transform 

Finally, updated approximation and detail coefficients are inverse Wavelet 
transformed to obtain the recovered signal. 

In this chapter the proposed denoising scheme was discussed in detail. This 
scheme was tested using three test signals. The test signal descriptions and the simulation 
results are presented in the next chapter. 
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IV. SIGNAL DESCRIPTION AND SIMULATION RESULTS 


The denoising scheme presented in the previous chapter was tested using three 
test signals. The codes used for the simulations are presented in Appendix A. In this 
chapter the test signals are described and the simulation results are summarized. 

A. SIGNAL DESCRIPTION 

Three different types of test signals are used: a sinusoid, a chirp with constant 
amphtude, and a chirp with amphtude increasing in an RC time constant fashion, not to 
exceed a given max im um value. These signals are described below. 

I. Sinusoidal Signal 

The frequency of the sinusoidal signal is selected randomly for every trial, where 
the frequency range is hmited to avoid abasing and DC signals. Figure 5 shows an 
example where /q = 10/512, and the samphng frequency f, = l. 



Figure 5. Sinusoidal test signal; frequency /g = 10/512, sampling frequency /^ = 1. 


Constant Amplitude Chirp 


The constant amphtude chirp used in the simulations is obtained by 

^ 3470 ^ 


s(t) - sin 


y40 + t j 


The signal sampled with samphng frequency /^ = 1 is shown in Figure 6.a. 


(4.1) 
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3. Chirp With an RC Time Constant-Like Amplitude Increase 

This signal has a frequency starting with a high value and decreasing in time. 
However, the amplitude of the signal starts with a small value and increases with time. 


The increasing amplitude chirp is obtained by 


s(t) = 


(t + 9)(5l2-t) . 271 547.05^ 

' -sin 


521 


t + 35.05 


(4.2) 


This test signal sampled with sampting frequency, f, = l, was obtained from Wavelab 
[13] and is illustrated in Figure 6.b. 




B. SIMULATION RESULTS 

MATLAB [14] simulations were performed to test the performance of the 
proposed schemes. One hundred trials are considered for both schemes for signal-to-noise 
ratio (SNR) values ranging between -6 dB and 6 dB. To perform the hypothesis test in 
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the second step of both FFT-based and WT-based schemes, the Bootstrap MATLAB 
Toolbox [15] was used. The MATLAB code is included in Appendix A. 


The mean square error (MSB) and cross-correlation coefficients between original 
and recovered signals were selected as performance criteria. The MSB was initially 
selected as it is commonly used in signal processing apphcations to measure signal 
differences in the time domain. It is defined as: 


1 'W f 1 V I ^ 2 


(4.3) 


where M is the number of trials, N is the signal length, Sj{i) and Sj{i) are the data 

sample of the original signal and the recovered signal attrial, respectively. Note that 
the MSB may not always be useful in evaluating actual performances. Bor example, in a 
denoising scheme when the recovered signal is close to the original version in most of the 
signal duration but it is very different for a short duration, the overall performance of the 
scheme can be satisfactory but at the same time the MSB may be large. In addition, 
simulations showed the MSB performances to be very similar on a large portion of the 
schemes investigated. Therefore, we considered an additional performance criterion, 
based on the cross-correlation coefficient to complement the information given by the 
MSB criterion. Recall that the norma liz ed cross-correlation coefficient is commonly used 
to evaluate signal similarities and is defined as: 

—r ^ 

V i=l i=l 


where () denotes complex conjugation. Note that the cross-correlation coefficient r 
should be equal to 1 for a perfectly reconshucted noise-free signal. The closer the 
magnitude of r gets to 0, the worse the denoised signal will be. The distance measure 
between 1 and the normalized cross-correlation coefficient is given by 


d 


J_ 

M 




j=i 


(4.5) 
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where M is the number of trials and r. is the cross-correlation at lag zero for the trial. 

The resulting measure d is the additional criterion used to evaluate the denoising 
schemes’ performances in extracting the noise-finee signal from the noisy signal. 

In our simulations sometimes the MSB performances of various schemes were 
close to each other and it was difficult to make a distinction among them. In this kind of 
cases the distance measure r was used in making a distinction. 

The proposed denoising schemes are compared with the original Wavelet 
shrinkage algorithm, which was introduced by Donoho and Johnstone [16, 17] to denoise 
signals embedded in additive white Gaussian noise with unit variance. Note that our 
schemes do not require knowledge of the noise variance levels to be apphed, which is not 
the case for the original Donoho and Johnstone scheme. Also note that variable noise 
variances had been used in the simulations for the proposed scheme. Therefore the setup 
for the simulations needed to be changed before comparing the simulations for the 
proposed schemes with the Wavelet shrinkage algorithm. The Wavelet shrinkage 
algorithm is defined and comparison results are presented in Appendix C. 

1. Fast Fourier Transform-Based Denoising 

Performance criteria for the three test signals are illustrated for both thresholding 
options, in Figures 7 to 9. Examples of noisy signals and their recovered versions 
obtained by FT-based denoising scheme are included in Appendix B. 
a. Sinusoidal Signal 

Results, shown in Figure 7, indicate that hard and soft thresholding 
schemes have similar performances at all SNR values considered. 
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(a) Distance measure for sinusoidal signal 



SNR 

Figure 7. Distance measure d and MSE performance criteria for the sinusoidal signal; FFT-based 

denoising. 


b. Constant Amplitude Chirp 

The MSE results shown in Figure 8 indicate no significant differences in 
thresholding scheme performances, while the distance measure results indicate better 
performance for the soft thresholding scheme at low and medium SNR levels. 


25 





(a) Distance measure for constant amplitude chirp signal 



SNR 

Figure 8. Distance measure d and MSB performance criteria for constant amplitude chirp signal; 

FBI-based denoising. 

c. Chirp With an RC Time Constant-Like Amplitude Increase 

Results shown in Figure 9 indicate that both thresholding schemes have 
similar MSB performances for SNR values below -3 dB, and that the soft thresholding 
scheme performs better for SNR values above -3 dB. The distance measure criterion 
indicates a better performance for soft thresholding for all SNR levels investigated. 
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(a) Distance measure for increasing amplitude chirp signal 



SNR 

Figure 9. Distance measure d and MSE performance criteria for increasing amplitude chirp signal; 

FFl-based denoising. 

2. Wavelet Transform-Based Denoising 

Many parameters affect the performance of the wavelet transform-based 
denoising scheme. Recall that we consider four thresholding implementation schemes on 
the approximation coefficients: approximation coefficients left unperturbed (method 1), 
approximation coefficients thresholded identical to the detail coefficients (method 2), 
approximation coefficients median filtered with a filter of length three (method 3), and 
using an order-two predictor (method 4). Methods 3 and 4 do not perform as well as the 
first two methods because they oversmooth the approximation coefficients sequence, 
causing loss of signal power and degradation in the performance of the denoising scheme. 

Other parameters are the number of decomposition levels and the order of the Daubechies 
wavelet to be used (3,4 and 5). The numbers in the figure legends next to the 
thresholding types are the Daubechies wavelet orders. MSE values obtained for aU 
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schemes are elose to eaeh other and do not provide mueh information about relative 
performanees. As a result, we eonsider the distanee measure d to eompare performanees. 
In addition, we only show results obtained for the best of the 4 methods for the 
approximation eoeffieients and the best of the 3 wavelet based deeompositions (3,4 or 5 
level deeompositions). Examples for noisy signals and their reeovered versions obtained 
with WT-based denoising sehemes are ineluded in Appendix B. 

a. Sinusoidal Signal 

Results show that best performanees are obtained with method 2 for this 
signal type. Figure 10 presents the WT-based seheme performance obtained for a 4-level 
deeomposition with different wavelet orders and both thresholding sehemes. Note that the 
sinusoidal signal has a eonstant frequeney remaining in one of the frequeney bands 
formed by the wavelet transform for the whole signal duration. Results show that the 
performanee of the WT-based denoising seheme depends on the frequeney of the speeifie 
test signal and the level of wavelet deeomposition. Note that test signal frequeneies were 
pieked randomly for eaeh trial. It should also be noted that when using method 1, the 
signal is left untouehed only when it is loeated in the lowest frequeney band (i.e., the 
band eontaining the approximation eoeffieients). However, noise when present is also left 
untouehed in that lowest frequeney band, while no sueh distinetion is present in method 
2 . 

Simulation results show that the soft thresholding seheme outperforms the 
hard thresholding implementation at SNR values below 2 dB, while aU thresholding 
sehemes perform similarly for higher SNRs. 
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(a) Distance measure for sinusoidal signal, 4 level decomposition, method 2 



SNR 

Figure 10. Distance measure and MSE for test signal 1, using a Wavelet-based 4-level 
decomposition and method 2, with Daubechies wavelet orders 3, 4 and 5. 


b. Constant Amplitude Chirp 

Results for the WT-based denoising scheme using 4-level decomposition 
and method 1 are presented in Figure 11. Method 1 was shown here because simulations 
indicated that method 1 has significantly higher performance than method 2 on the 
constant amphtude chirp signal. This result was to be expected as this signal has most of 
its power at low frequencies. Frequency components exist in several frequency bands, 
including the one containing the approximation coefficients. When method 2 is used, the 
sequence of approximation coefficients is tested for Gaussianity and may be estimated as 
Gaussian when the signal components are very short, resulting in thresholding of the 
signal components. However, such a problem does not exist in method 1 where the 
approximation coefficients are left untouched. Finally, simulation results show that a 
four-level decomposition performs better than a three or five-level decomposition. 
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(a) Distance measure for constant amplitude chirp signal, 4 level decomposition, method 1 



Figure 11. Performance for constant amplitude chirp signal, using a Wavelet based 4-level 
decomposition and method 1, with Daubechies wavelet orders 3, 4 and 5. 


c. Chirp With an RC Time Constant-Like Amplitude Increase 

Simulations show method 1 has the best performance of the 4 methods 
investigated and associated results are presented in Figure 12. Results shown in Figure 12 
indicate that similar performances are obtained for all thresholding schemes and wavelet 
orders considered. 


30 





(a) Distance measure for increasing ampiitude chirp signai, 4 ievei decomposition, method 1 



-6 -4 -2 0 2 4 6 

SNR 


Figure 12. Performance for increasing amplitude chirp signal, using a Wavelet-based 4-level 
decomposition and method 1, with Daubechies wavelet orders 3, 4 and 5. 


In this chapter the test signals were described and the simulation results 
were summarized. The conclusions are presented in the next chapter. 
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V. CONCLUSIONS AND FUTURE WORK 


A. CONCLUSIONS 

This work considered a Bootstrap-based denoising scheme to denoise 
deterministic signals embedded in additive Gaussian noise. The Bootstrap is a technique 
derived for improving the accuracy of a parameter estimator, used especially in situations 
where conventional techniques are not vahd since the data considered is short in length. 
This technique is used for parameter estimation and hypothesis testing in our work. 

The proposed Bootstrap-based denoising scheme has four steps. First, the noisy 
signal is transformed into a two-dimensional domain, using the Short-Time Fast Fourier 
or the Wavelet Transform. Next, Gaussianity tests are performed on the transformed data 
on a frequency band basis. No processing is applied to the data when the tested data is 
found to be non-Gaussian, as this indicates signal components in the data are dominant 
over Gaussian noise components. However, denoising is applied when the data tested is 
found to be Gaussian, which indicates that the noise components are dominant. Denoising 
is obtained by thresholding the frequency components in specific frequency bands to 
minimize noise effects. Finally, the recovered signal is obtained by applying the 
appropriate inverse transform. 

Denoising scheme performances were investigated using three test signals: a 
sinusoid, a constant amphtude chirp and a chirp with increasing amphtude. MSB and 
cross-correlation measures were selected to investigate the relative performances of all 
four denoising schemes considered in this work. In most cases, for a given SNR the MSB 
values obtained for schemes with different decomposition types and levels were similar at 
all SNRs and were not useful in discriminating between the various schemes investigated. 
However, the distance measure was more sensitive and showed more discriminating 
information. Results show that FFT-based schemes perform better than WT-based 
schemes on the stationary sinusoid signal type, whereas WT-based scheme outperforms 
the BBT-based scheme on the chirp signal types considered. In the FBI case, soft and 
hard thresholding schemes perform similarly on the sinusoid, while soft thresholding 
outperforms hard thresholding on the chirp signals. Results also show that soft 
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thresholding performs better on the stationary sinusoidal signal type but performs similar 
to hard thresholding on the chirps in the WT case. Note that thresholding methods are 
performed on the detail coefficients in the WT-based denoising case. However, using 
them directly on the approximation coefficients may cause signal loss. Therefore, four 
schemes were considered in the simulations: Leaving the approximation coefficients 
untouched, using the same thresholding schemes as the detail coefficients, using a median 
filter of order 3 on the approximation coefficients, and using a predictor of order 2 on the 
approximation coefficients. The second method performs best for the sinusoidal test 
signal whereas the first method performs the best for the chirp type test signals. 

Daubechies wavelets of orders 3, 4 and 5 were selected in the WT-based 
denoising scheme. Simulations show that wavelets of this family with higher orders did 
not improve performance for signals with the given length. 

The proposed denoising scheme was compared with the Wavelet shrinkage 
algorithm, which was introduced by Donoho and Johnstone [16, 17] to denoise signals 
embedded in additive white Gaussian noise with unit variance. Results indicate that the 
proposed WT-based denoising scheme performs better than the wavelet shrinkage 
algorithm for the sinusoidal test signal. However, the wavelet shrinkage algorithm 
outperforms the proposed scheme for the chirp type test signals. It should be noted that 
these results (Donoho and Johnstone) are restricted to the case of noise with unit 
variance. 

B. FUTURE WORK 

Our study was restricted to one wavelet family (Daubechies) only and further 
investigations should consider other types of wavelets to investigate the impact of the 
specific wavelet type on the resulting denoising performances. 


34 



APPENDIX A. MATLAB CODES 


MATLAB simulations were performed to were used to test the performance of 
the denoising scheme. The codes used are presented in this chapter. The first and second 
sections present the codes for the FFT- and WT-based denoising schemes, respectively. 
The third section presents the sub-functions that were used in both of these codes. 


A. FFT-BASED DENOISING 
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testkurtvV: 

implement s 
three test 

the FFT-based 
signals 

denoising scheme 

on % 
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o. 

SYNTAX : 

testkurtvV 



o 

o, 

o 

o. 

INPUT : 

none 



o 

o, 

o 

o. 

OUTPUT : 

simulation 

results saved 

on disk 

o 

o, 

o 

o. 

SUB FUNC : 

MakeSignal.m 
kurtosistest.m 


o 

o, 

o 

o, 

o 

o, 

o 


Written by Hasan E. KAN 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


clear 

clc 

close all 
t=l:512; 

ampl=sqrt(10."( [ .6 0 -.6])); 

kurlim=4.5; 

for sigtype=l:3 

for trial=l:100 

recksoft = zeros(length(ampl) , 512 ) ; 
reckhard=zeros(length(ampl) , 512) ; 
if sigtype==l 

s = sin(2*pi*t*(1 + ceil(rand*240))/512) ; 
s = s/sqrt(mean(s."2) ) ; 

sigtypes='sin' ; 
elseif sigtype==2 
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s=sin(34.7./[.01:.01:612]); 
s = s(41:552);s = s/sqrt(mean(s."2) ) ; 
sigtypes='cch'; 
elseif sigtype==3 

s=makesignal('Doppler' , 521) ; 
s = s(10:521);s = s/sqrt(mean(s."2) ) ; 

sigtypes='ech' ; 

end 

signal=ones(length(ampl),1)*s; 
sigpwr=sum(signal.^2,2)/512; 

n=randn(length(ampl),512); 

noise=((ampl./std(n'))'*ones(1,512)).*n; 

noipwr=sum(noise.^2,2) / 512; 

SNR=10*logl0(sigpwr./noipwr); 

x=signal+noise; 

for snr=l:length(ampl)%1:21 
fx=zeros(61,32); 
for segmentno=l:61 

fx(segmentno,:)=fft(x(snr,(segmentno- 
1)*8 + 1: (segmentno-1)*8 + 32) .*triang(32) ') ; 
end 

fksoft=zeros(61,32); 
fksoftr=zeros(61,17) ; 
fksofti=zeros(61,17) ; 
fkhard=zeros (61,32) ; 
fkhardr=zeros (61,17); 
fkhardi=zeros (61,17); 
kurr=zeros(length(ampl),17); 
kuri=zeros(length(ampl),17); 
for segmentno=l:17 

[fksoftr(:,segmentno),fkhardr(:,segmentno),kurr(segmentno)]=kurto 
sistest (real(fx(:,segmentno)),kurlim) ; 

if imag(fx(:,segmentno))~=0 

[fksofti(:, segmentno) , fkhardi(:,segmentno),kuri(segmentno)]=kurto 
sistest(imag(fx(:,segmentno)),kurlim); 
end 

end 

fksoft(:,l:17)=fksoftr(:,l:17)+j*fksofti ( :,1:17); 
fkhard(:,1:17)=fkhardr(:,1:17)+j*fkhardi(:, 1:17) ; 
fksoft(:,18:32)=conj(fliplr(fksoft(:,2:16))); 
fkhard(:,18:32)=conj(fliplr(fkhard(:, 2:16) ) ) ; 
for segmentno=l:61 

recksoft(snr, (segmentno-1)*8 + 1: (segmentno- 
1)*8+32)=recksoft(snr,(segmentno-1)*8+1:(segmentno- 
1)*8+32)+ifft(fksoft(segmentno,:)); 
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reckhard(snr,(segmentno-1)*8+1:(segmentno- 
1)*8+32)=reckhard(snr,(segmentno-1)*8+1:(segmentno- 
1)*8+32)+ifft(fkhard(segmentno,:)); 
end 

end 

recksoft=recksoft*0.5; 
reckhard=reckhard*0.5; 

msesoft(:,trial)=sum((signal-recksoft).^2,2)/512; 
msehard(:,trial)=sum((signal-reckhard)."2,2)/512; 
for snr=l:21 

corsoft(snr,trial)=xcorr(s,recksoft(snr,:),0,'coeff'); 

corhard(snr,trial)=xcorr(s,reckhard(snr, :),0, 'coeff') ; 
end 

end 

save(sprint f('mse%s',sigtypes), 'msemed', 'msesoft', 'msehard') 
save(sprintf('cor%s',sigtypes),'corsoft','corhard') 

end 


B. WT-BASED DENOISING 


2 - 2 - 2 - 2 - 2 - 2 - 2 - 9 - 2 - 2 - 2 - 9 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 9 - 9 - 2 - 9 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 9 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


testkurtvV: 

; implements the WT-based denoising scheme on 
three test signals 

SYNTAX : 

; testwav2 

INPUT : 

; none 

OUTPUT : 

; simulation results saved on disk 

SUB FUNC : 

: MakeSignal.m 
kurtosistest.m 

Written by 

Hasan E. KAN 


2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 9 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 9 - 2 - 9 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 9 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


clear 

clc 

close all 
t=l:512; 

ampl=sqrt(10.^( [ .6 0 -.6])); 
kurlim=4.5; 

rec5db3soft=zeros(length(ampl),512,4); 
rec5db3hard=zeros(length(ampl) ,512,4); 
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rec4db3soft = zeros(length(ampl) ,512,4); 
rec4db3hard=zeros(length(ampl) ,512,4); 
rec3db3soft = zeros(length(ampl) ,512,4); 
rec3db3hard=zeros(length(ampl) , 512,4) ; 

rec5db4soft = zeros(length(ampl) , 512, 4 ) ; 
rec5db4hard=zeros(length(ampl) , 512,4) ; 
rec4db4soft = zeros(length(ampl) , 512, 4 ) ; 
rec4db4hard=zeros(length(ampl) , 512,4 ) ; 
rec3db4soft = zeros(length(ampl) , 512, 4) ; 
rec3db4hard=zeros(length(ampl) ,512,4); 

rec5db5soft = zeros(length(ampl) ,512,4); 
rec5db5hard=zeros(length(ampl) ,512,4); 
rec4db5soft = zeros(length(ampl) ,512,4); 
rec4db5hard=zeros(length(ampl) ,512,4); 
rec3db5soft = zeros(length(ampl) ,512,4); 
rec3db5hard=zeros(length(ampl) , 512, 4) ; 

for sigtype=l:3 

for trial=l:100 
if sigtype==l 

s = sin(2*pi*t*(1 + ceil(rand*2 40))/512);s = s/sqrt(sum(s."2)/512) ; 
sigtypes='sin'; 
elseif sigtype==2 

s = sin(34.7./[.01: .01:612]);s = s(41:552);s = s/sqrt(sum(s.^2)/512) ; 
sigtypes='cch' ; 
elseif sigtype==3 

s=makesignal('Doppler',521) ;s = s(10:521) ;s = s/sqrt(sum(s."2)/512) ; 
sigtypes='ech' ; 

end 

signal=ones(length(ampl),1)*s; 
sigpwr=sum(signal.^2,2) / 512; 

n=randn(length(ampl),512); 

noise=((ampl./std(n'))'*ones(1,512)).*n; 

noipwr=sum(noise.^2,2) / 512; 

SNR=10*logl0(sigpwr./noipwr) ; 

x=signal+noise; 

for snr=l:length(ampl) 

[cdb3(snr,:),ldb3]=wavedec(x(snr,:),3,'db3'); 
a3db3(snr,:)=appcoef(cdb3(snr,:),ldb3,'db3',3); 

[dldb3(snr,:),d2db3(snr,:),d3db3(snr,:)]=detcoef(cdb3(snr,:),ldb3 
,[1,2,3]); 

[a4db3(snr,:),d4db3(snr,:)]=dwt(a3db3(snr,:),'db3'); 
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[a5db3(snr,:),d5db3(snr,:)]=dwt(a4db3(snr,:),'db3'); 


[cdb4(snr,:),ldb4]=wavedec(x(snr,:),3,'db4'); 
a3db4(snr, :)=appcoef(cdb4(snr, :),ldb4, 'db4',3); 

[dldb4(snr,:),d2db4(snr,:),d3db4(snr,:)]=detcoef(cdb4(snr,:),ldb4 
, [1,2,3] ) ; 

[a4db4(snr, :),d4db4 (snr, :)]=dwt(a3db4(snr, :), 'db4') ; 
[a5db4(snr,:),d5db4(snr,:)]=dwt(a4db4(snr,:),'db4'); 

[cdb5(snr, :),ldb5]=wavedec(x(snr, :),3,'db5'); 
a3db5(snr,:)=appcoef(cdb5(snr,:),ldb5,'db5',3); 

[dldbS(snr,:),d2db5(snr,:),d3db5(snr,:)]=detcoef(cdb5(snr,:),IdbS 
, [1,2,3] ) ; 

[a4db5(snr,:),d4db5(snr,:)]=dwt(a3db5(snr,:),'dbS'); 
[aSdbS(snr,:),d5db5(snr,:)]=dwt(a4db5(snr,:),'dbS'); 

end 

ra5db3soft=zeros(length(ampl),length(a5db3));ra5db3hard=zeros(len 
gth(ampl),length(a5db3)); 

ra4db3soft=zeros(length(ampl),length(a4db3)) ; ra4db3hard=zeros(len 
gth(ampl),length(a4db3) ) ; 

ra3db3soft=zeros(length(ampl),length(a3db3));ra3db3hard=zeros(len 
gth(ampl),length(a3db3)); 

rd5db3soft=zeros(length(ampl),length(d5db3(1,:)));rd5db3hard=zero 
s (length(ampl),length(d5db3(1, :)) ) ; 

rd4db3soft=zeros(length(ampl),length(d4db3));rd4db3hard=zeros(len 
gth(ampl),length(d4db3)); 

rd3db3soft=zeros(length(ampl),length(d3db3));rd3db3hard=zeros(len 
gth(ampl),length(d3db3)); 

rd2db3soft=zeros(length(ampl),length(d2db3));rd2db3hard=zeros(len 
gth(ampl),length(d2db3)); 

rdldb3soft=zeros(length(ampl),length(dldb3));rdldb3hard=zeros(len 
gth(ampl),length(dldb3) ) ; 


ra5db4soft=zeros(length(ampl),length(aSdb4));ra5db4hard=zeros(len 
gth(ampl),length(a5db4)); 

ra4db4soft=zeros(length(ampl),length(a4db4));ra4db4hard=zeros(len 
gth(ampl),length(a4db4) ) ; 

ra3db4soft = zeros(length(ampl),length(a3db4)) ; ra3db4hard=zeros(len 
gth(ampl),length(a3db4) ) ; 
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rd5db4soft=zeros(length(ampl),length(d5db4));rd5db4hard=zeros(len 
gth(ampl),length(d5db4)); 

rd4db4soft=zeros(length(ampl),length(d4db4));rd4db4hard=zeros(len 
gth(ampl),length(d4db4)); 

rd3db4soft = zeros(length(ampl),length(d3db4)) ; rd3db4hard=zeros(len 
gth(ampl),length(d3db4) ) ; 

rd2db4soft=zeros(length(ampl),length(d2db4));rd2db4hard=zeros(len 
gth(ampl),length(d2db4)); 

rdldb4soft=zeros(length(ampl),length(dldb4));rdldb4hard=zeros(len 
gth(ampl),length(dldb4)); 

ra5db5soft=zeros(length(ampl),length(aSdbS)) ; ra5db5hard=zeros(len 
gth(ampl),length(aSdbS) ) ; 

ra4db5soft=zeros(length(ampl),length(a4db5));ra4db5hard=zeros(len 
gth(ampl),length(a4db5)); 

ra3db5soft=zeros(length(ampl),length(a3db5));ra3db5hard=zeros(len 
gth(ampl),length(a3db5)) ; 

rd5db5soft=zeros(length(ampl),length(dSdbS));rd5db5hard=zeros(len 
gth(ampl),length(dSdbS)); 

rd4db5soft = zeros(length(ampl),length(d4db5)) ; rd4db5hard=zeros(len 
gth(ampl),length(d4db5) ) ; 

rd3db5soft=zeros(length(ampl),length(d3db5));rd3db5hard=zeros(len 
gth(ampl),length(d3db5)); 

rd2db5soft=zeros(length(ampl),length(d2db5));rd2db5hard=zeros(len 
gth(ampl),length(d2db5)); 

rdldb5soft=zeros(length(ampl),length(dldbS));rdldb5hard=zeros(len 
gth(ampl),length(dldbS) ) ; 
for method=l:4 

if method==l 

ra5db3soft=a5db3;ra5db3hard=a5db3; 
ra4db3soft=a4db3;ra4db3hard=a4db3; 
ra3db3soft=a3db3;ra3db3hard=a3db3; 
ra5db4 soft=a5db4;ra5db4hard=a5db4; 
ra4db4 soft=a4db4;ra4db4hard=a4db4; 
ra3db4 soft=a3db4;ra3db4hard=a3db4; 
raSdbSsoft=a5db5;ra5db5hard=a5db5; 
ra4db5soft=a4db5;ra4db5hard=a4db5; 
ra3db5soft=a3db5;ra3db5hard=a3db5; 
methodl='no operation on approximation coeff.'; 
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elseif method==2 

for snr=l:length(ampl) 

[raSdbSsoft(snr,:),raSdbShard(snr,:),kuraSdbS(snr)]=kurtosistest( 
a5db3(snr, :),kurlim) ; 

[ra4db3soft(snr,:),ra4db3hard(snr,:),kura4db3(snr)]=kurtosistest( 
a4db3(snr, :),kurlim) ; 

[ra3db3soft(snr, :),ra3db3hard(snr, :),kura3db3(snr)]=kurtosistest ( 
a3db3(snr,:),kurlim); 

[ra5db4soft(snr,:),ra5db4hard(snr,:),kura5db4(snr)]=kurtosistest( 
a5db4(snr, :),kurlim) ; 

[ra4db4soft(snr,:),ra4db4hard(snr,:),kura4db4(snr)]=kurtosistest( 
a4db4(snr, :),kurlim) ; 

[ra3db4soft(snr, :),ra3db4hard(snr, :),kura3db4(snr)]=kurtosistest ( 
a3db4(snr,:),kurlim); 

[raSdbSsoft(snr,:),raSdbShard(snr,:),kuraSdbS(snr)]=kurtosistest( 
aSdbS(snr, :),kurlim) ; 

[ra4db5soft(snr,:),ra4db5hard(snr,:),kura4db5(snr)]=kurtosistest( 
a4db5(snr, :),kurlim) ; 

[ra3db5soft(snr,:),ra3db5hard(snr,:),kura3db5(snr)]=kurtosistest( 
a3db5(snr,:),kurlim); 

end 

methodl='using threshold on approximation 

coeff. ' ; 

elseif method==3 

ra5db3soft=medfiltl(a5db3' , 3) ';ra5db3hard=ra5db3soft; 

ra4db3soft=medfiltl(a4db3',3) ';ra4db3hard=ra4db3soft; 

ra3db3soft=medfiltl(a3db3',3)';ra3db3hard=ra3db3soft; 

ra5db4soft=medfiltl(a5db4',3)';ra5db4hard=ra5db4soft; 

ra4db4soft=medfiltl(a4db4',3)';ra4db4hard=ra4db4soft; 

ra3db4soft=medfiltl(a3db4',3)';ra3db4hard=ra3db4soft; 

ra5db5soft=medfiltl(aSdbS',3)';ra5db5hard=ra5db5soft; 

ra4db5soft=medfiltl(a4db5',3)';ra4db5hard=ra4db5soft; 

ra3db5soft=medfiltl(a3db5',3)';ra3db5hard=ra3db5soft; 

methodl='using medfilt3 on approximation coeff.'; 
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else 


rra5db3=xcorr(a5db3',2,'unbiased');rra5db3=rra5db3(3:5,4*[0:2]+l) 

r 

rra4db3=xcorr(a4db3',2, 'unbiased');rra4db3=rra4db3(3:5,4*[0:2]+l) 

r 

rra3db3=xcorr(a3db3',2,'unbiased');rra3db3=rra3db3(3:5,4*[0:2]+!) 

r 

rra5db4=xcorr(a5db4',2,'unbiased');rra5db4=rra5db4(3:5,4*[0:2]+l) 

r 

rra4db4=xcorr(a4db4',2,'unbiased');rra4db4=rra4db4(3:5,4*[0:2]+l) 

r 

rra3db4=xcorr(a3db4',2,'unbiased');rra3db4=rra3db4(3:5,4*[0:2]+l) 

r 

rra5db5=xcorr(a5db5',2, 'unbiased');rra5db5=rra5db5(3:5, 4*[0:2]+!) 

r 

rra4db5=xcorr(a4db5',2, 'unbiased');rra4db5=rra4db5(3:5,4*[0:2]+l) 

r 

rra3db5=xcorr(a3db5',2,'unbiased');rra3db5=rra3db5(3:5,4*[0:2]+l) 

r 

for snr=l:length(ampl) 

Rxa5db3=toeplitz(rra5db3(:,snr));Rxxa5db3=Rxa5db3(1:2, 1:2) ; 

Rxa4db3=toeplitz(rra4db3(:,snr) );Rxxa4db3=Rxa4db3(1:2, 1:2) ; 

Rxa3db3=toeplitz(rra3db3(:,snr));Rxxa3db3=Rxa3db3(1:2,1:2); 

Rxa5db4=toeplitz(rra5db4(:,snr));Rxxa5db4=Rxa5db4(1:2,1:2); 

Rxa4db4=toeplitz(rra4db4(:,snr));Rxxa4db4=Rxa4db4(1:2,1:2); 

Rxa3db4=toeplitz(rra3db4(:,snr));Rxxa3db4=Rxa3db4(1:2,1:2); 

Rxa5db5=toeplitz(rra5db5(:,snr) );Rxxa5db5=Rxa5db5(1:2, 1:2) ; 

Rxa4db5=toeplitz(rra4db5(:,snr));Rxxa4db5=Rxa4db5(1:2,1:2); 

Rxa3db5=toeplitz(rra3db5(:,snr));Rxxa3db5=Rxa3db5(1:2,1:2); 

aa5db3=Rxxa5db3\Rxa5db3(2:3,1); 
aa4db3=Rxxa4db3\Rxa4db3(2:3,1); 
aa3db3=Rxxa3db3\Rxa3db3(2:3,1); 
aa5db4=Rxxa5db4\Rxa5db4(2:3,1); 
aa4db4=Rxxa4db4\Rxa4db4(2:3,1); 
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aa3db4=Rxxa3db4\Rxa3db4(2:3,1) ; 
aa5db5=Rxxa5db5\Rxa5db5(2:3,1) ; 
aa4db5=Rxxa4db5\Rxa4db5(2:3,1) ; 
aa3db5=Rxxa3db5\Rxa3db5(2:3,1) ; 





ra5db3soft 

(snr. 

: )=filter([ 

0; 

aa5db3] 

, 1, a 5 db 3 

(snr. 

:) ); 







ra4db3soft 

(snr. 

: )=filter([ 

0; 

aa4db3] 

, 1, a 4 db 3 

(snr. 

:) ); 







ra3db3soft 

(snr. 

: )=filter([ 

0; 

aa3db3] 

, 1, a 3 db 3 

(snr. 

:) ); 







ra5db4 soft 

(snr. 

: )=filter([ 

0; 

aa5db4] 

, 1, a 5 db 4 

(snr. 

:) ); 







ra4db4soft 

(snr. 

: )=filter([ 

0; 

aa4db4] 

, 1, a 4 db 4 

(snr. 

:) ); 







ra3db4 soft 

(snr. 

: )=filter([ 

0; 

aa3db4] 

, 1, a 3 db 4 

(snr. 

:) ); 







raSdbSsoft 

(snr. 

: )=filter([ 

0; 

aaSdbS] 

, 1, a 5 db 5 

(snr. 

:) ); 







ra4db5soft 

(snr. 

: )=filter([ 

0; 

aa4db5] 

, 1, a 4 db 5 

(snr. 

:) ); 







ra3db5soft 

(snr. 

: )=filter([ 

0; 

aa3db5] 

, 1, a 3 db 5 

(snr. 

:) ); 






end 







ra5db3soft(:, 1: 

: 2) =a5db3 ( :, 1:2) 

r 



ra4db3soft(:, 1: 

; 2) =a4db3 ( :, 1:2) 

r 


ra3db3soft(:,l:2)=a3db3(:,1:2) ; 
ra5db4soft(:,l:2)=a5db4(:,1:2) ; 
ra4db4soft(:,l:2)=a4db4(:,1:2) ; 
ra3db4soft(:,l:2)=a3db4(:,1:2) ; 
raSdbSsoft(:,l:2)=a5db5(:,1:2) ; 
ra4db5soft(:,l:2)=a4db5(:,l:2) ; 
ra3db5soft(:,l:2)=a3db5(:,l:2) ; 
ra5db3hard=ra5db3soft; 
ra4db3hard=ra4db3soft; 
ra3db3hard=ra3db3soft; 
ra5db4hard=ra5db4 soft; 
ra4db4hard=ra4db4 soft; 
ra3db4hard=ra3db4 soft; 
ra5db5hard=ra5db5soft; 
ra4db5hard=ra4db5soft; 
ra3db5hard=ra3db5soft; 

methodl='prediction on approximation coeff.'; 
end 

for snr=l:length(ampl) 

[rd5db3soft(snr,:),rd5db3hard(snr,:),kurd5db3(snr)]=kurtosistest( 
d5db3(snr, :),kurlim) ; 

[rd4db3soft(snr,:),rd4db3hard(snr,:),kurd4db3(snr)]=kurtosistest( 
d4db3(snr, :),kurlim) ; 
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[rd3db3 soft 

(snr,:),rd3db3hardi 

(snr. 


kurd3db3 

(snr)] 

=kurtosistest 

d3db3(snr,: 

),kurlim); 






[rd2db3 soft 

(snr,:),rd2db3hardi 

(snr. 

:) , 

kurd2db3 

(snr)] 

=kurtosistest 

d2db3(snr,: 

),kurlim); 






[rdldb3 soft 

(snr, :),rdldb3hardi 

(snr. 

:) , 

kurdldbS 

(snr)] 

=kurtosistest 

dldb3(snr,: 

),kurlim); 






[rd5db4 soft 

(snr,:),rd5db4hardi 

(snr. 

:) , 

kurdSdb4 

(snr)] 

=kurtosistest 

d5db4(snr,: 

),kurlim); 






[rd4db4 soft 

(snr,:),rd4db4hardi 

(snr. 

:) , 

kurd4db4 

(snr)] 

=kurtosistest 

d4db4(snr,: 

),kurlim); 






[rd3db4 soft 

(snr,:),rd3db4hardi 

(snr. 

:) , 

kurd3db4 

(snr)] 

=kurtosistest 

d3db4(snr,: 

),kurlim); 






[rd2db4 soft 

(snr,:),rd2db4hardi 

(snr. 

:) , 

kurd2db4 

(snr)] 

=kurtosistest 

d2db4(snr,: 

),kurlim); 






[rdldb4 soft 

(snr,:),rdldb4hardi 

(snr. 

:) , 

kurdldb4 

(snr)] 

=kurtosistest 

dldb4(snr,: 

),kurlim); 






[rdSdbS soft 

(snr,:),rdSdbShardi 

(snr. 

:) , 

kurdSdbS 

(snr)] 

=kurtosistest 

dSdbS(snr,: 

),kurlim); 






[rd4db5soft 

(snr,:),rd4db5hardi 

(snr. 

:) , 

kurd4dbS 

(snr)] 

=kurtosistest 

d4db5(snr,: 

),kurlim); 






[rd3db5 soft 

(snr, :),rd3dbShardi 

(snr. 

:) , 

kurdSdbS 

(snr)] 

=kurtosistest 

d3db5(snr,: 

),kurlim); 






[rd2db5 soft 

(snr,:),rd2dbShardi 

(snr. 

:) , 

kurd2dbS 

(snr)] 

=kurtosistest 

d2db5(snr,: 

),kurlim); 






[rdldbS soft 

(snr,:),rdldbShardi 

(snr. 


kurdldbS 

(snr)] 

=kurtosistest 

dldbS(snr,: 

),kurlim); 







end 







RSdb3soft=[raSdbSsoft 

rdSdbSsoft 

rd4db3soft 

rd3db3soft 

rd2db3soft rdldbSsoft]; 






RSdb3hard=[raSdbShard 

rdSdbShard 

rd4db3hard 

rd3db3hard 

rd2db3hard rdldbShard]; 






R4db3soft=[ra4db3soft 

rd4db3soft 

rd3db3soft 

rd2db3soft 

rdldbSsoft]; 







R4db3hard=[ra4db3hard 

rd4db3hard 

rd3db3hard 

rd2db3hard 

rdldbShard]; 







R3db3soft=[ra3db3soft 

rd3db3soft 

rd2db3soft 


rdldbSsoft 
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rdldb3hard] 

R3db3hard=[ra3db3hard 

f 

rd3db3hard 

rd2db3hard 

rd3db4 soft 

R5db4soft=[ra5db4 soft 
rd2db4soft rdldb4soft]; 

rd5db4 soft 

rd4db4 soft 

rd3db4hard 

R5db4hard=[ra5db4hard 
rd2db4hard rdldb4hard]; 

rd5db4hard 

rd4db4hard 

rd2db4 soft 

R4db4soft=[ra4db4 soft 
rdldb4 soft]; 

rd4db4 soft 

rd3db4 soft 

rd2db4hard 

R4db4hard=[ra4db4hard 
rdldb4hard]; 

rd4db4hard 

rd3db4hard 

rdldb4 soft] 

R3db4soft=[ra3db4 soft 

f 

rd3db4 soft 

rd2db4 soft 

rdldb4hard] 

R3db4hard=[ra3db4hard 

f 

rd3db4hard 

rd2db4hard 

rd3db5soft 

R5db5soft=[ra5db5soft 
rd2db5soft rdldb5soft]; 

rd5db5soft 

rd4db5soft 

rd3db5hard 

R5db5hard=[ra5db5hard 
rd2db5hard rdldb5hard]; 

rd5db5hard 

rd4db5hard 

rd2db5soft 

R4db5soft=[ra4db5soft 
rdldb5soft]; 

rd4db5soft 

rd3db5soft 

rd2db5hard 

R4db5hard=[ra4db5hard 
rdldb5hard]; 

rd4db5hard 

rd3db5hard 

rdldb5soft] 

R3db5soft=[ra3db5soft 

r 

rd3db5soft 

rd2db5soft 


R3db5hard=[ra3db5hard 

rd3db5hard 

rd2db5hard 


rdldbShard]; 

for snr=l:length(ampl) 

rec5db3soft(snr, method)=waverec(R5db3soft(snr, :),[ 20 20 36 68 
131 258 512],'db3'); 

rec5db3hard(snr,:,method)=waverec(R5db3hard(snr,:),[20 20 

36 68 131 258 512],'db3'); 

rec4db3soft(snr,:,method)=waverec(R4db3soft(snr,:),[36 36 68 131 
258 512],'db3'); 

rec4db3hard(snr, :,method)=waverec(R4db3hard(snr,:), [36 36 

68 131 258 512],'db3'); 

rec3db3soft(snr, :,method)=waverec(R3db3soft(snr, :),[ 68 68 131 258 
512],'db3'); 

rec3db3hard(snr,:,method)=waverec(R3db3hard(snr,:),[68 68 

131 258 512],'db3'); 


rec5db4soft(snr,:,method)=waverec(R5db4soft(snr,:),[22 22 38 70 
133 259 512],'db4'); 
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rec5db4hard(snr,:,method)=waverec(R5db4hard(snr,:),[22 22 

38 70 133 259 512] , 'db4'); 

rec4db4soft(snrmethod)=waverec(R4db4soft(snr38 38 70 133 
259 512], 'db4') ; 

rec4db4hard(snr, :,method)=waverec(R4db4hard(snr,:), [38 38 

70 133 259 512] , 'db4'); 

rec3db4soft(snr,:,method)=waverec(R3db4soft(snr,:),[70 70 133 259 
512],'db4'); 

rec3db4hard(snr,:,method)=waverec(R3db4hard(snr,:),[70 70 

133 259 512],'db4'); 


rec5db5soft(snr,:,method)=waverec(R5db5soft(snr,:),[24 24 40 71 
134 260 512] , 'db5' ) ; 

rec5db5hard(snr, :,method)=waverec(R5db5hard(snr,:), [24 24 

40 71 134 260 512] , 'db5'); 

rec4db5soft(snr,:,method)=waverec(R4db5soft(snr,:),[40 40 71 134 
260 512], 'db5') ; 

rec4db5hard(snr, :,method)=waverec(R4db5hard(snr,:), [40 40 

71 134 260 512] , 'db5'); 

rec3db5soft(snr, :,method)=waverec(R3db5soft(snr, :),[71 71 134 260 
512],'db5'); 

rec3db5hard(snr,:,method)=waverec(R3db5hard(snr,:),[71 71 

134 260 512] , 'db5') ; 

end 


end 

signa=repmat(signal, [1,1,4]); 


mse5db3soft(:,trial,:)=sum( 
rec5db3soft).^2,2)/512; 

mse5db3hard(:,trial,:)=sum( 
rec5db3hard).^2,2)/512; 

mse4db3soft(:,trial,:)=sum( 
rec4db3soft).^2,2)/512; 

mse4db3hard(:,trial,:)=sum( 
rec4db3hard)."2,2)/512; 

mse3db3soft(:,trial,:)=sum( 
rec3db3soft).^2,2)/512; 

mse3db3hard(:,trial,:)=sum( 
rec3db3hard)."2,2)/512; 


(signa- 
(signa- 
(signa- 
(signa- 
(signa- 
(signa- 
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mse5db4soft(:,trial,:)=sum((signa- 
rec5db4soft).^2,2)/512; 

mse5db4hard(:,trial,:)=sum((signa- 
rec5db4hard).^2,2)/512; 

mse4db4soft(:,trial,:)=sum((signa- 
rec4db4soft).^2,2)/512; 

mse4db4hard(:,trial,:)=sum((signa- 
rec4db4hard)."2,2)/512; 

mse3db4soft(:,trial,:)=sum((signa- 
rec3db4soft).^2,2)/512; 

mse3db4hard(:,trial,:)=sum((signa- 
rec3db4hard)."2,2)/512; 

mse5db5soft(:,trial,:)=sum((signa- 
recSdbSsoft).^2,2)/512; 

mseSdbShard(:,trial,:)=sum((signa- 
recSdbShard).^2,2)/512; 

mse4db5soft(:,trial,:)=sum((signa- 
rec4db5soft).^2,2)/512; 

mse4db5hard(:,trial,:)=sum((signa- 
rec4db5hard)."2,2)/512; 

mse3db5soft(:,trial,:)=sum((signa- 
rec3db5soft)."2,2)/512; 

mse3db5hard(:,trial,:)=sum((signa- 
rec3db5hard)."2,2)/512; 
for snr=l:21 

for method=l:4 

cor3db3soft(snr,trial,method)=xcorr(s,rec3db3soft(snr, :, method) , 0 
,'coeff'); 

cor3db4soft(snr,trial,method)=xcorr(s,rec3db4soft(snr,:,method),0 
,'coeff'); 

cor3db5soft(snr,trial,method)=xcorr(s,rec3db5soft(snr, :,method),0 
,'coeff'); 

cor4db3soft(snr,trial,method)=xcorr(s,rec4db3soft(snr,:,method),0 
,'coeff'); 

cor4db4soft(snr,trial,method)=xcorr(s,rec4db4soft(snr,:,method),0 
,'coeff'); 

cor4db5soft(snr,trial,method)=xcorr(s,rec4db5soft(snr, :, method) , 0 
,'coeff'); 

cor5db3soft(snr,trial,method)=xcorr(s,rec5db3soft(snr,:,method),0 
,'coeff'); 

cor5db4soft(snr,trial,method)=xcorr(s,rec5db4soft(snr, :,method),0 
,'coeff'); 
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corSdbSsoft(snr,trial,method)=xcorr(s,rec5db5soft(snr, method) , 0 
,'coeff'); 

cor3db3hard(snr,trial,method)=xcorr(s,rec3db3hard(snr,:,method),0 
,'coeff'); 

cor3db4hard(snr,trial,method)=xcorr(s, rec3db4hard(snr, :,method),0 
,'coeff'); 

cor3db5hard(snr,trial,method)=xcorr(s,rec3db5hard(snr,:,method),0 
,'coeff'); 

cor4db3hard(snr,trial,method)=xcorr(s,rec4db3hard(snr,:,method),0 
,'coeff'); 

cor4db4hard(snr,trial,method)=xcorr(s,rec4db4hard(snr, : , method) , 0 
,'coeff'); 

cor4db5hard(snr,trial,method)=xcorr(s,rec4db5hard(snr,:,method),0 
,'coeff'); 

cor5db3hard(snr,trial,method)=xcorr(s, rec5db3hard(snr, :,method),0 
,'coeff'); 

cor5db4hard(snr,trial,method)=xcorr(s,rec5db4hard(snr,:,method),0 
,'coeff'); 

corSdbShard(snr,trial,method)=xcorr(s,recSdbShard(snr,:,method),0 
,'coeff'); 

end 

end 

end 


save(sprintf('mse%s31vl',sigtypes) 
3hard','mse3db4soft','mse3db4hard' 

save(sprintf('mse%s41vl',sigtypes) 
3hard','mse4db4soft','mse4db4hard' 


'msemed','mse3db3soft','mse3db 
'mse3db5soft','mse3db5hard') 

'msemed','mse4db3soft','mse4db 
'mse4db5soft','mse4db5hard') 


save(sprintf('mse%s51vl',sigtypes),'msemed','mse5db3soft','mseSdb 
3hard','mse5db4soft','mse5db4hard','mseSdbSsoft','mseSdbShard') 

save(sprintf('cor%s31vl',sigtypes),'cor3db3soft','cor3db3hard','c 
or3db4soft','cor3db4hard','cor3db5soft','cor3db5hard') 


save(sprintf('cor%s41vl',sigtypes),'cor4db3soft','cor4db3hard','c 
or4db4soft','cor4db4hard','cor4db5soft','cor4db5hard') 

save(sprintf('cor%s51vl',sigtypes),'cor5db3soft','cor5db3hard','c 

or5db4soft','cor5db4hard','corSdbSsoft','corSdbShard') 

end 
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C. SUB FUNCTIONS 
1. MakeSignal 

This code is part of Wavelab version 802 [11]. 


function sig = MakeSignal(Name,n) 

% MakeSignal — Make artificial signal 
% Usage 

% sig = MakeSignal(Name,n) 

% Inputs 

% Name string: 'HeaviSine', 'Bumps', 'Blocks', 

% 'Doppler', 'Ramp', 'Cusp', 'Sing', 'HiSine', 

% 'LoSine', 'LinChirp', 'TwoChirp', 'QuadChirp', 

% 'MishMash', 'WernerSorrows' (Heisenberg), 

% 'Leopold' (Kronecker), 'Piece-Regular' (Piece-Wise 


Smooth), 

% 'Riemann','HypChirps','LinChirps', 'Chirps', 'Gabor' 

% 'sineoneoverx','Cusp2','SmoothCusp','Gaussian' 

% 'Piece-Polynomial' (Piece-Wise 3rd degree polynomial) 

% n desired signal length 

% Outputs 

% sig 1-d signal 

o, 

o 

% References 

% Various articles of D.L. Donoho and I.M. Johnstone 

o, 

o 


if nargin > 1, 
t = (1:n) ./n; 
end 

if strcmp(Name,'HeaviSine'), 
sig = 4.*sin (4*pi.*t); 

sig = sig - sign(t - .3) - sign(.72 - t); 
elseif strcmp(Name,'Bumps'), 


pos = 

[ .1 

.13 .15 

.23 

.25 

.40 

.44 

. 65 

.76 

.78 

.81] ; 

hgt = 

[ 4 

5 3 

4 

5 4 . 

,2 2. 

, 1 4 . 

, 3 

3.1 5, 

, 1 4 . 

.2] ; 

wth = 

[ .005 

.005 . 

006 

.01 . 

,01 , 

,03 , 

, 01 

.01 . 

,005 

.008 


.005] ; 

sig = zeros (size (t)); 
for j =1:length(pos) 

sig = sig + hgt(j)./( 1 + abs((t - 
pos ( j)) ./wth(j) )) .^4; 
end 

elseif strcmp(Name,'Blocks'), 

pos = [ .1 .13 .15 .23 .25 .40 .44 .65 .76 .78 .81]; 

hgt = [4 (-5) 3 (-4) 5 (-4.2) 2.1 4.3 (-3.1) 2.1 (- 

4.2) ] ; 

sig = zeros (size (t)); 
for j=l:length(pos) 

sig = sig + (1 + sign(t-pos ( j))).*(hgt ( j)/2) ; 

end 
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elseif strcmp(Name,'Doppler'), 

sig = sqrt(t.*(1-t) ) .*sin ( (2*pi*l.05) ./(t+.05)); 

elseif strcmp(Name,'Ramp'), 
sig = t - (t >= .37); 
elseif strcmp(Name,'Cusp'); 

sig = sqrt(abs(t - .37)); 
elseif strcmp(Name,'Sing'), 
k = floor(n * . 31 ); 
sig = 1 ./abs(t - (k+.5)/n); 
elseif strcmp(Name,'HiSine'), 

sig = sin( pi * (n * .6902) .* t); 

elseif strcmp(Name,'LoSine'), 

sig = sin( pi * (n * .3333) .* t); 

elseif strcmp(Name,'LinChirp'), 

sig = sin(pi .* t .* ( (n .* .500) .* t)); 

elseif strcmp(Name,'TwoChirp'), 

sig = sin (pi .* t .* (n .* t)) + sin((pi/3) .* t .* (n 

.* t) ) ; 

elseif strcmp(Name,'QuadChirp'), 

sig = sin( (pi/3) .* t .* (n .* t.^2)); 

elseif strcmp(Name,'MishMash'), % QuadChirp + LinChirp + 

HiSine 


sig = 

sin 

( (pi/3) .* 

t . * 

(n . * t.' 

^2) ) ; 

sig = 

sig 

+ sin ( pi 

* (n 

* .6902) 

.* t) ; 

sig = 

sig 

+ sin (pi , 

. * t . 

,* (n .* 

. 125 . * t) ) ; 


elseif strcmp(Name,'WernerSorrows'), 

sig = sin( pi .* t .* (n/2 .* t."2)) ; 

sig = sig + sin( pi * (n * .6902) .* t); 

sig = sig + sin (pi .* t .* (n .* t)); 

pos = [ .1 .13 .15 .23 .25 .40 .44 .65 .76 .78 .81]; 

hgt =[4 5 3 4 5 4.22.14.3 3.15.14.2]; 

wth = [.005 .005 .006 .01 .01 .03 .01 .01 .005 .008 

.005] ; 

for j =1:length(pos) 

sig = sig + hgt (j)./( 1 + abs((t - 
pos (j)) ./wth( j) )) .^4; 
end 

elseif strcmp(Name,'Leopold'), 

sig = (t == floor(.37 * n)/n); % Kronecker 

elseif strcmp(Name,'Riemann'), 
sqn = round(sqrt(n)); 

sig = t .* 0; % Riemann's Non-differentiable Function 

sig( (1:sqn) .^2) = 1. ./ (l:sqn); 

sig = real (ifft(sig)); 

elseif strcmp(Name,'HypChirps'), % Hyperbolic Chirps of 
Mallat's book 

alpha= 15*n*pi/1024; 
beta = 5*n*pi/1024; 

t = (1.001:l:n+.001)./n; 
f1 = zeros (1,n); 

f2 = zeros (1,n); 

fl = sin(alpha./(.8-t)).*(0.l<t).*(t<0.68) ; 
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f2 = sin (beta./(.8-t)) .* (0 . l<t) . * (t<0.75); 

M = round(0.65*n); 

P = floor(M/ 4 ); 

enveloppe = ones(l,M); % the rising cutoff function 
enveloppe (1 : P) = (1 + sin(-pi/2+((1:P)-ones (1,P)) . / (P- 

1)*pi))/2; 

enveloppe(M-P + 1:M) = reverse(enveloppe(1:P)) ; 
env = zeros(1, n) ; 

env(ceil(n/10):M+ceil(n/10)-1) = enveloppe(1:M); 
sig = (fl+f2).*env; 

elseif strcmp(Name,'LinChirps'), % Linear Chirps of 
Mallat's book 

b = 100*n*pi/1024; 

a = 250*n*pi/1024; 

t = (1 :n) . /n; 

A1 = sqrt( (t-l/n) .* (1-t)); 

sig =Al.*(cos((a*(t).^2)) +cos((b*t+a*(t)."2))); 
elseif strcmp(Name,'Chirps'), % Mixture of Chirps of 
Mallat's book 

t = (l:n)./n.*10.*pi; 

fl = cos(t."2*n/1024); 

a = 30*n/1024; 

t = (1:n)./n.*pi; 

f2 = cos(a.* (t. "3)); 

f2 = reverse (f2); 

ix = (-n:n)./n.*20; 

g = exp(-ix."2*4*n/1024); 

11 = (n/2+1:n/2+n); 

12 = (n/8+1:n/8+n); 

j = (l:n)/n; 

f3 = g(11).*cos(50.*pi.*j*n/1024); 
f4 = g(12).*cos(350.*pi.*j*n/1024); 
sig = fl+f2+f3+f4; 

enveloppe = ones(l,n); % the rising cutoff function 
enveloppe(1:n/8) = (1 + sin (-pi/2+((1:n/8)- 
ones(l,n/8))./(n/8-l)*pi))/2; 

enveloppe(7*n/8+l:n) = reverse(enveloppe(1:n/8)); 
sig = sig.*enveloppe; 

elseif strcmp(Name,'Gabor'), % two modulated Gabor 
functions in 


% Mallat's book 


N = 512; 
t = (-N:N)*5/N; 
j = (1:N)./N; 
g = exp(-t.^2*20); 

11 = (2*N/4+l:2*N/4+N); 

12 = (N/4+1:N/4+N); 

sigl = 3*g(il).*exp(i*N/l6.*pi.*j); 
sig2 = 3*g(12).*exp(i*N/4.*pi.*j); 
sig = sigl+sig2; 

elseif strcmp(Name,'sineoneoverx'), % sin(l/x) in Mallat's 


book 
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N = 1024; 
i = (-N+1:N); 
i(N) = 1/100; 
i = i./(N-1); 
sig = sin(1.5./(i)) ; 
sig = sig(513:1536); 
elseif strcmp(Name,'Cusp2'), 

N = 64; 
i = (1:N) ./N; 

X = (l-sqrt(i)) + i/2 -.5; 

M = 8*N; 

sig = zeros(1,M); 
sig(M-1.5.*N+1:M-.5*N) = x; 
sig(M-2.5*N+2:M-1.5.*N+1) = reverse(x); 
sig(3*N+1:3*N + N) = .5*ones(1,N); 
elseif strcmp(Name,'SmoothCusp'), 
sig = MakeSignal('Cusp2'); 

N = 64; 

M = 8*N; 
t = (1:M)/M; 
sigma = 0.01; 
g = exp (-. 5 . * (abs (t- 
.5)./sigma)."2)./sigma./sqrt(2*pi); 
g = fftshift(g); 
sig2 = iconv(g',sig)'/M; 
elseif strcmp(Name,'Piece-Regular'), 
sigl=-l5*MakeSignal('Bumps',n); 
t = (1:fix(n/12)) ./fix(n/12); 

sig2=-exp(4*t); 
t = (l:fix(n/7)) ./fix(n/7); 

sig5=exp(4*t)-exp(4); 
t = (l:fix(n/3)) ./fix(n/3); 

sigma=6/40; 

sig6=-70*exp(-( (t-1/2) .*(t-l/2) )/(2*sigma^2)); 
sig(l:fix(n/7))= sig6(l:fix(n/7)); 

sig((fix(n/7)+l):fix(n/5))=0.5*sig6((fix(n/7)+l):fix(n/5)); 

sig( (fix(n/5)+l) :fix(n/3))=sig6 ( (fix(n/5)+l) :fix(n/3)); 

sig( (fix(n/3)+l) :fix(n/2))=sigl ( (fix(n/3)+l) :fix(n/2)); 
sig((fix(n/2)+l):(fix(n/2)+fix(n/12)))=sig2; 
sig((fix(n/2)+2*fix(n/12)):- 
1: (fix(n/2)+fix(n/12)+l))=sig2; 

sig(fix(n/2)+2*fix(n/12)+fix(n/20)+1:(fix(n/2)+2*fix(n/12)+3*fix( 
n/20)))=... 

-ones(1,fix(n/2)+2*fix(n/12)+3*fix(n/20)-fix(n/2)-2*fix (n/12)- 
fix(n/20))*25; 

k=fix(n/2)+2*fix(n/12)+3*fix(n/20); 
sig( (k+1) : (k+fix(n/7)))=sig5; 
diff=n-5*fix(n/5); 

sig(5*fix(n/5)+1:n)=sig(diff:-l:1); 
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% zero-mean 
bias=sum(sig)/n; 
sig=bias-sig; 

elseif strcmp(Name,'Piece-Polynomial'), 
t = (l:fix(n/5)) ./fix(n/5); 

sigl=20*(t."3+t.^2+4); 
sig3=40*(2.*t."3+t) + 100; 
sig2=10.*t.^3 + 45; 
sig4=16*t."2+8.*t+16; 
sig5=20*(t+4); 

sig6(l:fix(n/10))=ones(l,fix(n/10)); 
sig6=sig6*20; 
sig(l:fix(n/5))=sigl; 

sig(2*fix(n/5):-l:(fix(n/5)+l))=sig2; 

sig( (2*fix(n/5)+l) :3*fix(n/5))=sig3; 

sig((3*fix(n/5)+l):4*fix(n/5))=sig4; 

sig( (4*fix(n/5)+1) :5*fix(n/5))=sig5(fix(n/5) :-l: 1); 

diff=n-5*fix(n/5); 

sig(5*fix(n/5)+1:n)=sig(diff:-l:1); 

%sig( (fix(n/20)+l) : (fix(n/20)+fix(n/10) ))=- 
ones(l,fix(n/10))*20; 

sig((fix(n/20)+l):(fix(n/20)+fix(n/10)))=ones(l,fix(n/10))* 

10; 

sig((n-fix(n/10)+l):(n+fix(n/20)- 
fix(n/10)))=ones (l,fix(n/20))*150; 

% zero-mean 
bias=sum(sig)/n; 
sig=sig-bias; 

elseif strcmp(Name,'Gaussian'), 
sig=GWN(n,beta); 
g=zeros(1,n); 
lim=alpha*n; 
mult=pi/(2*alpha*n); 
g(1:lim) = (cos(mult*(1:lim))) ."2 ; 
g( (n/2 + 1) :n)=g( (n/2) :-l:l); 
g = rnshift(g,n/2); 
g=g/norm(g); 
sig=iconv(g,sig); 

else 

disp(sprintf('MakeSignal: I don*t recognize 
<<%s>>',Name)) 

disp('Allowable Names are:') 
disp('HeaviSine'), 
disp('Bumps'), 
disp('Blocks'), 
disp ( 'Doppler'), 
disp('Ramp'), 
disp('Cusp'), 
disp('Crease'), 
disp('Sing'), 
disp('HiSine'), 
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o\® o\® o\® 


disp 

'LoSine'), 

disp 

'LinChirp'), 

disp 

'TwoChirp'), 

disp 

'QuadChirp'), 

disp 

'MishMash'), 

disp 

'WernerSorrows'), 

disp 

'Leopold'), 

disp 

'Sing'), 

disp 

'HiSine'), 

disp 

'LoSine'), 

disp 

'LinChirp'), 

disp 

'TwoChirp'), 

disp 

'QuadChirp'), 

disp 

'MishMash'), 

disp 

'WernerSorrows'), 

disp 

'Leopold'), 

disp 

'Riemann'), 

disp 

'HypChirps'), 

disp 

'LinChirps'), 

disp 

'Chirps ' ), 

disp 

'sineoneoverx'), 

disp 

'Cusp2 ' ), 

disp 

'SmoothCusp'), 

disp 

'Gabor'), 

disp 

'Piece-Regular'); 

disp 

'Piece-Polynomial'); 

disp 

'Gaussian'); 


end 


Originally made by David L. Donoho. 
Function has been enhanced. 


o, 

o 

% Part of WaveLab Version 802 
% Built Sunday, October 3, 1999 8:52:27 AM 
% This is Copyrighted Material 
% For Copying permissions see COPYING.m 
% Comments? e-mail wavelab@stat.stanford.edu 

o, 

o 
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2 , 


Kurtosistest 


kurtosistest 


implements thebootstrap-based hypothesis 
test H: kurtosis of the input=3 against 
kurtosis of input~=3 and performs soft and 
hard thresholding on the data 


SYNTAX : [soft,hard,kur] = kurtosistest[data,kurlim] 

INPUT : data: data sequence to be tested 

kurlim: kurtosis limit not to be exceeded 
thresholding 


OUTPUT 


soft: soft thresholded data ° 
hard: hard thresholded data ^ 
kur: kurtosis value used in soft thresholding^ 


SUB FUNC 


none 


Written by Hasan E. KAN 


- 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 9 - 2 - 2 - 2 - 2 - 9 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 9 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooo 


function [varargout]=kurtosistest(varargin); 

datain=varargin{1}; 
kurlim=varargin{2}; 

hardout=zeros(size(datain) ) ; 

test=boottest(datain, 'kurtosis',3,1, .05,99,25) ; 
if test==l 

softout=datain; 
hardout=datain; 
kur=-l; 

else 

storage=[];temp=datain; 

kur=mean(bootstrp(99,'kurtosis',datain)); 
while kur>kurlim 

if length(storage)>round(length(temp)/3) 
kur=kurlim; 

else 

locations=find(abs(temp)==max(abs(temp))) ' ; 
storage=[storage; locations]; 
temp (locations)=0; 

kur=mean(bootstrp(9 9, 'kurtosis',temp)) ; 

end 


end 
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softout=datain*abs(3-kur)/(kurlim-3); 
if -isempty(storage) 

softout(storage)=datain(storage); 
hardout(storage)=datain(storage); 

end 

end 

varargout{l} = softout; 
varargout{2} = hardout; 
varargout{3} = kur; 
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APPENDIX B. SIMULATION RESULTS 


Simulation results for the three test signal types are presented in detail in this 
Appendix. Descriptions for these test signals were given in Chapter 4. Results include 
performance criteria obtained for the MSB and distance measure defined in Chapter 4 for 
SNR levels between -6 to 6 dB, for one hundred trials per SNR level. Finally, the 
original, noisy and the recovered versions of one trial are shown for SNR values -6 dB, 

0 dB and 6 dB to provide a visual representation of the FFT-based denoising scheme. 

A. FFT-BASED DENOISING 
I. Sinusoidal Signal 

The frequency for the sinusoidal signal was selected randomly for each of the 100 
trials so as to generahze the simulation results. The frequency was lim ited to avoid 
aliasing or DC signals. Minimum and maximum allowed frequencies for the sinusoidal 
test signal in the simulations were 2/512 Hz and 241/512 Hz respectively, where the 
sampling frequency was selected as 1 Hz. Figure 13 shows the performance curves for 
FFT-based denoising scheme. These curves indicate that the thresholding schemes 
perform similarly on this type of signal. Figures 14 to 16 illustrate noisy and recovered 
versions of a sinusoidal signal at SNR values of -6 dB, 0 dB and 6 dB respectively. 
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Figure 13. FFT-based denoising for sinusoidal signal type. 












Original signal 



Figure 14. FFT-based denoising for sinusoidal signal type, SNR=-6 dB. 
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Original signal 



Figure 15. FFT-based denoising for sinusoidal signal type, SNR=0 dB. 
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Original signal 



Figure 16. FFT-based denoising for sinusoidal signal type, SNR=6 dB. 
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2. Constant Amplitude Chirp Signal 

Soft and hard thresholding perform similarly according to MSB. However, the 
distance measure shows that soft thresholding outperforms hard thresholding, especially 
at low and medium SNRs. This result can be seen in Figure 17. Figures 18 to 20 illustrate 
noisy and recovered versions of a sinusoidal signal at SNR values of -6 dB, 0 dB and 
6 dB respectively. 


(a) Distance measure for constant amplitude chirp signal 



(b) MSE for constant amplitude chirp signal 



SNR 

Figure 17. FFT-based denoising for constant amplitude chirp type signal. 
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Original signal 



Figure 18. FFT-based denoising for constant amplitude chirp type, SNR=-6 dB 






Original signal 


2 

0 
-2 
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Noisy signal (SNR=0dB) 
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0 
-2 
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Recovered signal with soft thresholding 
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0 
-2 

0 50 100 150 200 250 , 300, 350 400 450 500 

Recovered signal with hard thresnolding 
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0 
-2 

0 50 100 150 200 250 300 350 400 450 500 






Figure 19. FFT-based denoising for constant amplitude chirp type, SNR=0 dB. 
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Original signal 
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0 
-2 

0 50 100 150 200 250 , 300, 350 400 450 500 

Recovered signal with hard thresnolding 




Figure 20. FFT-based denoising for constant amplitude chirp type, SNR=6 dB. 
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3. Chirp with an Amplitude that Increases in an RC Time Constant 
Fashion 

The results shown in Figure 21 indicate that the thresholding schemes have 
similar MSE performances at SNR values lower than -3 dB, whereas soft thresholding 
performs better at higher SNRs. According to the distance measure d, soft thresholding is 
better at all SNR levels. Figures 22 to 24 illustrate noisy and recovered versions of a 
sinusoidal signal at SNR values of -6 dB, 0 dB and 6 dB respectively. 


(a) Distance measure for increasing ampiitude chirp signal 




SNR 

Figure 21. FFT-based denoising for increasing amplitude chirp type signal. 
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Original signal 



0 50 100 150 200 250 300, 350 400 450 500 

Recovered signal with hard thresnolding 



Figure 22. FFl-based denoising for increasing amplitude chirp type, SNR=-6 dB. 
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Original signal 



0 50 100 150 200 250 , 300, 350 400 450 500 

Recovered signal with hard thresnolding 



Figure 23. FFT-based denoising for increasing amplitude chirp type, SNR=0 dB. 


Original signal 



Figure 24. FFT-based denoising for increasing amplitude chirp type, SNR=6 dB. 


69 




B. WT-BASED DENOISING 

Many different parameters affect the performance of the WT-based denoising 
scheme such as the number of decomposition levels, processing method for the 
approximation coefficients and the order of the Daubechies wavelet used. Results show 
no significant difference among the Daubechies wavelet orders considered. Only best 
results for each signal are illustrated in this section. Results shown include MSE and 
distance measure performances for SNR levels between -6 to 6 dB, where one hundred 
trials are used for each SNR level. Noisy and recovered versions of one trial are shown 
for SNR values -6 dB, 0 dB and 6 dB to provide a visual representation of the WT-based 
denoising scheme. 

1. Sinusoidal Signal 

Soft thresholding outperforms hard thresholding at low and medium SNR levels 
as shown in Figure 25. However, it does not perform as good as any of the thresholding 
methods used in the FFT-based denoising case for this signal type. Noisy and recovered 
versions of a sinusoidal signal at SNR values of -6 dB, 0 dB and 6 dB are shown in 
Figures 26 to 28. 
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(a) Distance measure for sinusoidal signal, 4 level decomposition, method 2 



-6 -4 -2 0 2 4 6 


(b) MSE for sinusoidal signal, 4 level decomposition, method 2 



-6 -4 -2 0 2 4 6 

SNR 


Figure 25. WT-based denoising for sinusoidal signal type, 4 level, method 2, with Daubechies 

wavelet orders 3, 4 and 5. 
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Figure 26. 


Noisy signal (SNR=-6dB) 



WT-based de noising for sinusoidal signal type, 4 level, method 2, Daubechies wavelet 4, 

SNR=-6 dB. 
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2 . 


Constant Amplitude Chirp 


Results in Figure 29 show that soft and hard thresholding performances are 
similar for this signal when using WT-based denoising scheme. The WT-based scheme 
performs better than the FFT-based scheme on this signal. Figures 30 to 32 illustrate 
noisy and recovered versions of a sinusoidal signal at SNR values of -6 dB, 0 dB and 
6 dB respectively. 

(a) Distance measure for constant amplitude chirp signal, 4 level decomposition, method 1 




Figure 29. WT-based denoising for constant amplitude chirp type, 4 level, method 1, with 

Daubechies wavefet orders 3, 4 and 5. 
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Figure 31. WT-based denoising for constant amplitude chirp type, 4 level, method 1, Daubechies 

wavelet 4, SNR=0 dB. 


77 





























0 50 100 150 200 250 300 350 400 450 500 


Figure 32. WT-based denoising for constant amplitude chirp type, 4 level, method 1, Daubechies 

wavelet 4, SNR=6 dB. 
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3. Chirp with an Amplitude that Increases in an RC Time Constant 
Fashion 

Soft and hard thresholding methods perform similarly according to the 
performance criteria shown in Figure 33. WT-based denoising scheme outperforms FFT- 
based denoising scheme on this signal. Figures 34 to 36 illustrate noisy and recovered 
versions of a sinusoidal signal at SNR values of -6 dB, 0 dB and 6 dB respectively. 

(a) Distance measure for increasing amplitude chirp signal, 4 level decomposition, method 1 



-6 -4 -2 0 2 4 6 

SNR 


Figure 33. WT-based denoising for increasing amplitude chirp type, 4 level, method 1, with 

Daubechies wavelet orders 3, 4 and 5. 
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Original signal 



Figure 34. WT-based denoising for increasing amplitude chirp type, 4 level, method 1, Daubechies 

wavelet 5, SNR=-6 dB. 
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Figure 36. WT-based denoising for increasing amplitude chirp type, 4 level, method 1, Daubechies 

wavelet 5, SNR=6 dB. 
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APPENDIX C. COMPARISON WITH WAVELET SHRINKAGE 

ALGORITHM 


The proposed denoising schemes are compared with the original Wavelet 
shrinkage algorithm in this chapter. Wavelet shrinkage algorithm is a well-known 
denoising scheme. 


A. WAVELET SHRINKAGE ALGORITHM 

The Wavelet shrinkage algorithm was introduced by Donoho and Johnstone [16, 
17] to denoise signals embedded in additive white Gaussian noise with unit variance. 

This algorithm has three steps: 

a. Wavelet transformation. 


b. Noisy coefficients suppression by applying a non-linear thresholding 
technique. 

c. Inverse wavelet transformation. 

1. Choosing a Threshold Value 

The following four threshold selection schemes are available: 

a. Stein’s Unbiased Risk Estimator (SURE) Threshold 

The SURE threshold value is derived adaptively for each decomposition 
level by minimizing the Stein’s Unbiased Estimate of risk [18] for threshold estimates. 

b. Sqtwolog Threshold 

This method uses a fixed threshold value defined as: 

T = yj2log(length(signal)) . 

c. Heursure Threshold 


This method is a mixture of the preceding two methods. It uses the 
Sqtwolog threshold at low SNR levels and the SURE threshold at medium to high SNR 
levels. 
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d. Minimax Threshold 

This method uses a fixed threshold that gives minimax performance for 
the MSB. The minimax principle is used in statistics in order to design estimators, which 
realize the minimum of the maximum mean square error obtained for the worst function 
in a given set. So, minimax threshold is the value that provides the minimum MSB 
among the worst threshold values for the detail and approximation coefficient sets. 

2. Thresholding Methods 

a. Hard Thresholding 

This method sets the coefficients with absolute values below the chosen 
threshold to zero. 

b. Soft Thresholding 

This method first sets the coefficients with absolute values below a chosen 
threshold to zero. Then, it shrinks the remaining coefficients using the relationship 

c - sign(c)[|c|- T], where c is the coefficient to be thresholded, c is the thresholded 

coefficient, and T is the chosen threshold value. 

B. COMPARISONS 

Simulations were performed, using the wavelet shrinkage algorithm defined 
above, on the same signals that were used for the proposed scheme. However, different 
SNR values were obtained by using constant noise power and different values for the 
signal power in this experiment, contrary to what was done in the simulations considered 
for the proposed scheme in the main body of the thesis. The Heursure method was used 
to compute the threshold values for the wavelet shrinkage algorithm. A new simulation 
was performed with the proposed scheme under the same SNR conditions as well. 

Results shown in Bigures 37 to 39 indicate that the proposed WT-based denoising scheme 
performs better than the wavelet shrinkage algorithm for the sinusoidal test signal. 

However, the wavelet shrinkage algorithm outperforms the proposed scheme for the chirp 
type test signals. It should be noted that these results are restricted to the case of noise 
with unit variance. 
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1. Sinusoidal Signal 

Results shown in Figure 37 indicate that the MSB performances of the two 
denoising algorithms are similar, whereas the distance measure d shows that the proposed 
denoising scheme outperforms the wavelet shrinkage algorithm at low SNRs. 


(a) Distance measure for sinusoidal signal, 3 level decomposition 



SNR 

Figure 37. Wavelet shrinkage algorithm and proposed WT-based denoising scheme performances 
for the sinusoidal signal type, using 3-level decomposition and Daubechies wavelet 5. 
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2 . 


Constant Amplitude Chirp 


Both performance criteria show that the wavelet shrinkage algorithm outperforms 
the proposed denoising scheme for this signal. This can be seen in Figure 38. 


(a) Distance measure for constant amplitude chirp, 4 level decomposition 




SNR 

Figure 38. Wavelet shrinkage algorithm and proposed WT-based denoising scheme performances 
for the constant amplitude chirp type, using 4-level decomposition and Daubechies wavelet 4. 
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3. Chirp With an RC Time Constant-Like Amplitude Increase 

Both performance criteria show that the wavelet shrinkage algorithm outperforms 
the proposed denoising scheme for this signal as shown in Figure 39. 


(a) Distance measure for increasing amplitude chirp, 4 level decomposition 




SNR 

Figure 39. Wavelet shrinkage algorithm and proposed WT-based denoising scheme performances 
for the increasing amplitude chirp type, using 4-level decomposition and Daubechies wavelet 5. 


It should be noted that the new proposed method (based on kurtosis and Bootstrap 
method) seem to be inferior to the methods described in this appendix. But the methods 
in this appendix need additional information: the variance of the noise is unity, which 
may be an unreasonable constraint in some practical applications. 
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